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In  this  report  we  summarize  our  research  activities  during  the  three 
tnonth  period  February  i - April  30,  1975.  Our  research  program  started 
October  1,  1973*  During  the  first  year  we  concentrated  our  efforts  on  basic 
research  in  image  structure  analysis  and  image  modeling.  During  the  second 
year  we  are  continuing  our  basic  research,  but  in  the  meantime  also  look  into 
several  applications  in  the  information  extraction  area.  We  hope  that  on  th^ 
one  hand  the  results  of  our  basic  research  could  find  practical  applications, 
and  on  the  other  hand  problems  arising  from  practical  situations  would  guide  the 
direction  of  our  basic  research. 

In  basic  research,  our  emphasis  has  been  shifted  toward  image  understand- 
ing and  Information  extraction.  The  ultimate  goal  in  automatic  image  under- 
standing and  Information  extraction  is  to  build  a mechine  which  will  look  at  a 
scene  and  then  give  us  a qualitative  as  well  as  quantitative  description  of  it. 
This  formidable  task  can  be  broken  down  into  several  sub-tasks  as  indicated  by 
the  block  diagram  in  Fig.  1. 

The  sensor  looks  at  the  scene  and  outputs  image  data  which  are  generally 
multlspectral  and  multi  temporal . The  preprocessing  box  attempts  either  to  put 
the  image  into  a form  which  is  more  suitable  for  analysis  by  the  latter  boxes 
or  to  compress  the  image  data  for  transmission  to  a remote  location  or  storage 
for  future  use.  The  image  segmentation  box,  in  simple  instances,  locates 
objects  in  the  image  (for  example.  In  character  recognition,  it  locates  the 
characters).  For  more  complex  scenes.  It  segments  the  Image  into  regions. 

Each  of  these  regions  Is  classified  by  the  classification  box.  The  classifi- 
cation can  be  done  either  by  classical  decision-theoretic  methods  (feature 
extraction  followed  by  statistical  classification)  or  by  the  more  recent 
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Fig.  I Image  understanding  i.nd  Information  extraction 


syntactical  methods.  In  linguistic  terminology,  the  regions  are  primitives, 
and  the  classifier  finds  attributes  for  these  primitives.  Finally,  the 
structural  analysis  box  attempts  to  find  the  relationship  (spatial,  temporal, 
spectral)  among  the  regions  (primitives).  In  some  cases,  one  may  want  to 

develop  a grammar  for  the  primitives. 

The  functions  of  the  boxes  in  Fig.  1 overlap,  and  their  demarcation  is 
not  clear  cut.  Furthermore,  there  is  strong  feedback  in  the  system  (as  in- 
dicated by  the  dashed  arrows  in  the  block  diagram).  For  example,  structural 
analysis  may  reveal  that  certain  regions  are  classified  wrongly  and  thus 
force  a reclassification.  Similarly,  difficulties  encountered  In  classifi- 
cation may  make  it  necessary  to  reexamine  image  segmentation. 

The  research  projects  we  are  carrying  and  will  carry  out  fall  into  six 
overlapping  categories  corresponding  to  the  five  boxes  in  Fig.  1 and  applica- 
tions. in  the  present  report  our  emphasis  is  on  applications,  image  segmen- 
tation, and  pattern  classification.  We  mention  some  highlights  below. 

We  have  attacked  quite  successfully  four  application-oriented  problems, 
locating  large  airports  in  ERTS  imagery  (Huang  and  Chan),  detecting  water 
boundaries  in  ERTS  imagery  for  flood  maps  (Chen  and  Wintz),  change  detection 
in  missile  movie  frames  (Huang  and  Burnett),  and  TV  bandwidth  compression  for 

remotely  piloted  vehicle  images  (Wintz). 

In  image  segmentation  we  are  pursuing  two  approaches;  edge  detection, 
and  region  analysis.  In  the  case  of  edge  detection  we  are  particularly  in- 
terested in  developing  techniques  which  would  work  even  if  the  image  is  noisy 
and  smeared  (Tang  and  Huang,  Burnett  and  Huang).  The  investigation  on  the  use 
of  syntactic  approach  (Kelly,  in  the  last  Quarterly  Progress  Report)  is  being 

continued. 
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In  the  case  of  region  analysis  we  have  developed  earlier  the  BLOB 
algorithm  (WIntz  and  Gupta,  Semi-annual  Report.  Oct.  I,  1973  - March  31,  I974). 
We  are  currently  working  on  Improved  versions  of  BLOB,  and  In  the  meantime 
looking  Into  techniques  which  would  give  results  that  are  Independent  of  the 
order  of  processing.  One  such  technique  Is  the  valley-seeking  unsupervised 
clustering  method  of  Fukuncga  and  Narendra  developed  recently  at  Purdue.  The 
use  of  unsupervised  clustering  methods  for  Image  segmentation  looks  promising 
(Yoo  and  Huang).  The  success  of  region  analysis  depends  critically  on  a suit- 
able choice  of  feature  meaaureaenta.  The  most  Important  class  of  features  for 
this  application  Is  perhaps  texture  descriptors.  We  have  been  studying 
texture  descriptors  Intensively  and  have  developed  several  new  ones  which  give 
good  classification  results  (Mitchell  and  Myers). 

In  pattern  classification  one  Is  often  overwhelmed  by  a large  number  of 
potentially  useful  features.  On  the  other  hand,  one  can  usually  afford  to  use 
only  a few.  The  key  question  Is  how  to  select  m features  from  a potential  set 
of  n (>m)  In  an  optimal  way.  We  have  developed  a very  efficient  method  of 

feature  selection  which  Is  orders  of  magnitude  faster  than  exhaustive  search 
(Fukunaga  and  Narendra). 

Most  of  the  research  projects  reported  In  the  following  pages  are 
continuing. 
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LC''.ATE  AIRPORTS  IN  ERTS  IMAGERY 
T.  S.  Huang  and  C.  K.  Chan 


I . OBJECTIVE 

The  goal  of  this  project  is  to  develop  computer  techniques  for  locating 
large  airports  in  ERTS  mul t i spectral  imagery. 

1 1 . PROGRESS; 

Several  large  airports  in  the  United  States  (O'hare  Airport  in  Chicago, 

San  Francisco  Airport,  etc.)  are  picked  out  from  EATS  data  to  provide  a pre- 
liminary stud,  of  their  texture  and  properties.  It  is  found  that  in  all  four 
spectral  channels  used  by  ERTS.  the  spectral  reflectivity  of  the  airports  and 
highways  are  very  similar.  Thus  a preliminary  objective  is  to  pick  out  ard 
separate  the  airports  and  highways  (i.e.,  concrete)  from  the  other  contents  of 
the  ERTS  data.  A statistical  pattern  recognition  approach  is  used  in  this 
study  making  use  of  the  facilities  in  LARS.  A training  set  of  data  for  con- 
crete is  hand-picked  and  the  rest  of  thd  picture  is  classified  as  concrete  or 
not-concrete  from  the  statistics  of  the  training  set.  Fig.  1 shows  the  picture 
of  O'hare  Airport  in  channel  2 (wavelength  0.69  um)  and  Fig.  2 shows 
fled  and  thresholded  four  channels  output.  A threshold  of  t means  that  all 
points  having  a probability  of  less  than  t/100  of  being  concrete  is  thresholded 
away  (black  part  of  picture)  and  points  with  probability  > F/lOO  of  being  con 

Crete  shows  up  in  white  in  the  picture. 

Since  airports  and  highways  usually  are  composed  of  straight  lines,  a 

Fourier  transform  approach  is  used  to  process  the  classified  data.  A rourier 
transform  of  Fig.  2b  is  shown  in  Fig.  3.  The  Fourier  transform  is  then  filter- 
ed so  that  only  the  wedge  corresponding  to  the  parallel  runways  are  left. 
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O'hare  Airport 


Thresholded  inverse 
transform  of  Fig.  3 

Fig.  k 


Fourier  transform 
of  Fig.  2b 


Thresholded  single 
channel  (Ch.2)  output 


Inverse  transform  is  then  performed  and  the  results  thresholded.  See 
Fig.  A.  This  study  provides  an  Idea  of  how  well  one  can  Identify  airport 
runways  from  a picture  like  Fig.  I. 

The  Fourier  transform  approach  Is  also  applied  to  the  data  In  a single 
channel  (ch.  2).  The  thresholded  result  Is  shown  In  Fig.  5.  Although  It 
came  out  better  than  the  thresholded  k channels  classification  results,  we 

feel  that  a better  choice  of  training  set  would  greatly  Improve  the  result  In 
Fig.  I». 

III.  PROPOSAL  FOR  FURTHER  WORK 

So  far,  the  study  has  emphasized  removing  unwanted  Information  from  a 
picture  so  as  to  br.'ng  out  the  runways.  The  filter  used  for  the  Fourier  trans- 
form and  the  threshold  used  on  the  resulting  picture  are  handpicked.  Further 
work  will  be  aimed  at  automating  these  choices. 

The  first  thing  that  will  be  done  Is  to  threshold  the  Fourier  transform 
In  Fig.  3 so  that  only  the  prominent  lines  are  left.  The  Inverse  transform 
can  then  be  taken  so  that  hopefully  only  the  straight  lines  that  signify  the 
highways  and  runways  will  remain.  An  automatic  thresholding  scheme  can  then 
be  used  to  produce  results  similar  to  that  of  Fig.  5.  If  this  proves  suc- 
cessful, larger  sections  of  different  airports  will  be  tried. 
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BOUNDARY  DETECTION  FOR  U.S.  ARMY  CORPS  OF  ENGINEERS 

FLOOD  MAPS 

P.  H.  Ch«n  and  P.  A.  Wlntz 

The  U.S.  Army  Corps  of  Engineers  (USACE)  Is  In  charge  of  flood  control  In 
the  United  States.  In  order  to  determine  which  flood  gates  to  open  or  close 
and  In  order  to  predict  river  flood  conditions  they  use  mathematical  models. 
The  Inputs  to  the  mathematical  models  are  measurements  taken  frcm  river  water 
level  water  gauges,  rain  gauges,  snow  gauges,  etc.  There  Is  reason  to  believe 
that  the  data  required  by  their  models  could  be  obtained  In  a more  cost 
effective  manner  from  satellite  pictures.  For  example,  precipitation  maps  and 
maps  showing  changes  In  snow  cover  are  Indicative  of  the  amount  of  water  In- 
putted to  the  water  shed.  Furthermore,  flood  area  extent  maps  provide  re- 
quired data. 

The  purpose  of  this  project  Is  to  determine  flood  area  extent  maps  from 
ERTS  multlspectral  scanner  data  and  to  code  these  maps  In  an  efficient  manner 
In  order  to  reduce  the  bulk  of  the  data  Inputted  to  USACE  mathematical  models. 
A few  preliminary  results  are  Illustrated  In  Figs.  1-5.  Fig.  I shows  band  7 
of  an  ERTS  picture  encompassing  part  of  the  Ohio  River.  Water  shows  up  dark 
In  band  7 since  Infrared  energy  Is  absorbed  by  water.  Maximum  likelihood 
classifier  was  trained  on  some  river  and  non-river  data  and  then  used  to 
classify  each  picture  element  Into  either  water  or  not  water.  The  resulting 
classification  map  Is  presented  In  Fig.  2 where  water  was  made  bright  and  non- 
water picture  elements  dark.  Fig.  3 Illustrates  a clatalflcatlon  map  of  the 
same  area  when  there  Is  flood.  Figs,  k and  5 Illustrate  the  classification 
results  of  another  section  of  the  river  taken  at  two  different  times  (no 
flooding  and  flooding).  Preliminary  results  Indicate  that  2-dlmenslonal 
runlength  codes  and  contour  codes  will  significantly  reduce  the  amount  of  data 
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Fig.  3 Classification  Result  (Flooding) 
of  the  same  area  as  Flq.  I 


Fig.  k Classification  Result  (No  Flood) 
of  one  section  of  the  river 


tnat  must  be  processed  by  the  USAGE  mathematical  models.  This  is  extremely 
Important  since  the  Mississippi  Valley  watershed  covers  about  129  ERTS 
frames  consisting  of  approximately  30  million  bits. 


[1]  T.  S.  Huang,  W.  F.  Schreiber,  and  0.  J.  Tretiak,  "Contour  Coding  of 
''"«9es,  ^cture  Bandwidth  Compression,  ed.  by  T.  S.  Huang  and 
0*  J-  Tretiak,  pp.  1972. 
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EXPERIMENTS  IN  CHANGE  DETECTION  IN  MISSILE  IMAGERY 


T.  S.  Huang  and  J.  W.  Burnett 

I.  INTRODUCTION 

In  Figures  i (a)  and  (b)  we  have  two  successive  frames  (digitized)  of  a 
movie  film  of  a missile.  The  original  film  was  supplied  by  the  U.  S.  Army 
Vftilte  Sand  Missile  Range.  We  report  here  some  experiments  In  automatic 
detection  of  change  between  the  two  frames.  Change  detection  has  obvious 
applications  In  the  tracking  of  moving  objects. 

II.  CORRELATION  AND  DIFFERENCING 

Each  frame  was  digitized  to  128x128  samples,  and  8 bits/sample.  The  two 
frames  were  lined  up  by  cross*corre1at1on.  Then,  In  the  first  experiment, 
one  frame  was  subtracted  from  the  other  point  by  point,  resulting  In  Figure  2. 
Thresholding  Figure  2 yielded  Figure  3>  The  missile  motion  was  detected,  but 
the  result  contained  a considerable  amount  of  extraneous  points. 

III.  EDGE  COMPARISON 

In  the  second  experiment,  we  first  extract  the  edge  points  (using  a 
smoothed  gradient  operator)  In  the  two  frames,  yielding  Figures  4(a)  and  (b). 
Then  an  "exclusive-or"  operation  was  performed  (point  by  point)  on  these  two 
edge  pictures,  resulting  In  Figure  A point  In  Figure  5 was  white  If  and 
only  If  the  corresponding  points  In  Figs.  4(a)  and  (b)  were  of  different  colors 
(one  black  and  the  other  white).  Again  the  missile  motion  was  detected,  but 
again  we  have  a large  amount  of  extraneous  points. 

IV.  ELASTIC  EDGE  COMPARISON 

To  reduce  extraneous  points  In  the  result,  we  performed  a third  experi- 
ment In  which  the  two  edge  pictures  (Figures  4(a)  and  (b))  were  compared  with 
some  tolerance  allowed.  Specifically,  for  each  black  point  In  Figure  4(b), 
we  examined  the  3x3  point  neighborhood  around  the  corresponding  point  In 
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Figure  4(a).  We  painted  the  corresponding  point  in  Figure  6 black  if  and 
only  if  the  3x3  neighborhood  contained  no  black  points.  In  the  resulting 
Figure  6 we  observe  that  the  missile  mo*. ton  was  detected,  and  the  extraneous 
points  appeared  In  small  groups  (3  points  or  less  In  each  connected  group). 
These  extraneous  points  can  be  eliminated  completely  If  we  erase  all 
connected  point  groups  which  contain  3 or  less  points. 

V.  CONCLUSION 

The  elastic  edge  comparison  technique  seems  to  work  well  in  detecting 
changes  between  two  scenes  which  are  essentially  Identical  except  for  the 


changing  parts. 
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TV  BANDWIDTH  COMPRESSION  FOR  RPV's 
Paul  A.  Wintr. 

A critical  problem  In  the  RPV  system  Is  the  bandwidth  required  to  trans- 
mit the  TV  signals  from  the  vehicle  to  the  vehicle  controller.  There  are 
three  possibilities  for  reducing  the  bandwidth:  (I)  reduce  the  temporal 

resolution  by  transmitting  less  than  30  frames  per  second;  (2)  reduce  the 
spatial  resolution  by  coding  fewer  than  all  of  the  picture  elements  In  the 
fratries  transmitted;  (3)  reduce  the  spatial  redundancy  in  the  remaining  picture 
e I emen  ts . 

The  purpose  of  this  project  Is  to  Investigate  the  tradeoffs  between  the 
latter  two  methods  for  reducing  the  number  of  bits  required  to  code  the 
picture. 

• The  original  picture  was  digitized  Into  a 512x512  array  of  picture 
elements  of  8 bits  each.  The  resolution  was  then  reduced  to  256x256  picture 
elements  by  averaging  subsets  of  2x2  picture  elements  In  the  512x512  picture. 
The  512x512  picture  Is  presented  In  Fig.  I and  the  256x256  picture  Is  presented 
In  Fig.  5. 

Each  of  these  pictures  was  processed  by  the  hybrid  coding  scheme  due  to 
Habibl.  Each  line  of  picture  elements  was  first  divided  Into  sub-blocks 
containing  32  picture  elements  each.  Each  of  these  1x32  arrays  of  picture 
elements  was  transformed  with  the  cosine  transformation  to  obtain  32  co- 
efficients. This  cosine  transformation  was  performed  on  each  1x32  subsets  of 
the  entire  original  picture  In  order  to  reduce  the  redundancy  In  the  hori- 
zontal direction.  Successive  scan  lines  are  also  highly  correlated.  Since 
the  picture  elements  In  a 1x32  sub-array  of  picture  elements  are  highly 
correlated  with  the  picture  elements  in  the  1x32  sub-array  on  the  next  scan 
I me,  and  since  we  do  the  same  transformation  on  both  sub-arrays,  we  also 
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expect  the  coefficients  produced  by  the  transformations  to  be  highly  cor- 
related. This  redundancy  can  be  reduced  by  differencing  the  coefficients. 

This  can  be  accomplished  by  the  coding  scheme  known  a»  DPCM.  In  summary, 
the  cosine  transform  in  the  horizontal  direction  followed  by  DPCM  in  the 
vertical  direction  to  obtain  1x32  arrays  of  the  coefficient  differences. 

Since  the  coefficients  have  different  variances  it  follows  that  the 
differences  between  coefficients  will  also  have  different  variances.  In 
order  to  efficiently  code  the  coefficient  differences  we  use  block  quantization. 
That  Is,  we  use  a different  number  of  bits  to  code  the  different  coefficient 
differences.  The  bit  assignment  is  optimized  in  order  to  minimize  the  RMS 
quantization  error.  The  bit  allocations  are  listed  In  Table  1. 

The  picture  in  Fig.  1 was  coded  by  this  method  and  then  decoded  to  form 
a 512x512  replica  of  the  original  image.  The  decoded  pictures  for  2,  1,  and 
1/2  bits  per  picture  element  are  presented  in  Figs.  2,3»  end  4.  The  2-blt 
picture  has  only  a very  small  degradation  relative  to  the  original.  The  one- 
bit  picture  has  more  degradation.  The  1/2-blt  picture  not  only  has  very  notice- 
able degradation  but  the  edges  of  the  1x32  blocks  are  also  apparent. 

The  same  coding  scheme  was  applied  to  the  256x256  picture  of  Fig.  5.  The 
coded  pictures  were  decoded  to  form  a replica  of  the  256x256  image  which  was 
then  converted  to  a 512x512  image  by  repeating  each  picture  element.  The 
results  are  shown  in  Figs.  6,  7,  and  8 for  2,  1,  and  1/2  bits  per  picture 
element.  The  same  general  comments  apply  as  for  the  512x512  picture  except 
that  the  distortions  are  significantly  more  severe.  The  reason  for  this  Is 
that  there  is  less  correlation  between  the  256x256  picture  elements  than 
between  the  512x512  picture  elements  and  so  the  coding  algorithm  is  expected 
to  be  less  efficient. 

Finally,  we  compare  the  512x512  picture  coded  at  1/2-blt  per  picture 
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element  Illustrated  In  Fig.  k with  the  256x256  picture  coded  at  2 bits  per 
picture  element  shown  In  Fig.  8.  Both  of  these  pictures  require  the  same 
number  of  bits  to  transmit.  The  512x512  picture  at  1/2  bit  is  subjectively 
better  than  the  256x256  picture  at  2 bits. 

The  same  experiment  was  performed  on  a different  picture  and  the  decoded 
pictures  are  presented  In  Figs.  9 to  16.  Fig.  9 Is  the  512x512  original, 
and  Figs.  10,  II,  and  12  are  coded  at  2 bits,  I bit,  and  1/2  bit  respectively. 
Fig.  13  Is  the  256x256  original  and  Figs.  1^,  15,  and  16  correspond  to  coding 
this  picture  with  2 bits,  I bit,  and  1/2  bit.  Careful  observation  of  these 
pictures  leads  to  the  same  conclusion  as  for  the  previous  pictures. 

Reference 

A.  Habibl,  "Hybrid  Coding  of  Pictorial  Data,"  IEEE  Trans,  on  Communications, 
Vol.  COM-22,  No.  5,  PP.  6H-624,  May  197^. 
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Fig.  8 256x256  Coded  Image,  1/2  bit/picture  element 
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Fig.  |l»  256x256  Coded  Image,  2 bIts/pIcture  element 
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EDGE  EXTRACTION  TECHNIQUES 


G.  Tang  and  T.  S.  Huang 

I.  INTRODUCTION 

The  ability  of  a local  operator  to  extract  adequate  edge  information  is 
important  in  many  applications  of  image  processing.  It  has  been  used  In  image 
enhancement  [i]  and  image  segmentation. 

Various  local  edge  operators  have  been  proposed  [3] • [^] » [6],  [7]  since 
the  advent  of  computers.  Most  of  them  are  designed  originally  for  extracting 
edges  from  noise-free  pictures  or  nearly  no^se-free  pictures.  Based  on  the 
philosophy  behind  these  operators,  we  may  classify  them  into  2 categories; 
i.e.,  template-matching  operators  or  difference  operator.  Template  matching 
operators  report  the  closeness  between  an  ideal  edge  paradigm  and  an  actual 
local  region  on  the  picture  being  processed.  Difference  operators  report  the 
difference  in  brightness. 

A general  scheme  to  extr>'ict  edges  from  a noisy  degraded  picture  is  shown 

% 

in  Fig.  1.  The  pre-processor  is  used  to  reduce  the  Influence  of  random  dis- 
turbances to  a less  significant  degree  so  that  the  difference  operator  can 
report  meaningful  edge  Information.  Di fference  operator  reports  the  gray 
level  difference  in  a certain  way.  The  tester  is  used  to  decide  whether  the 
iterative  application  of  the  difference  operator  is  required  and  the  number  of 
iterations.  [7l  The  post-processor  has  two  functions.  The  first  is  to  emphasize 
the  edge  Information  reported  from  the  difference  operator  so  that  subtle 
edges  wi  i i not  be  missed  and  to  reduce  the  thickness  of  reported  edges.  The 
second  is  to  change  a mui tiple-ievei  picture  to  a binary  picture  by  setting  a 
threshold. 

In  this  report  we  proposed  an  edge-emphasizing  post-processor  which  is 
proven  to  improve  the  quality  of  the  edges  a great  deal.  We  are  going  to  call 
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this  processor  "Local  unit  transformer"  which  means  the  transform  of  gray 
levels  from  a global  unit  to  a local  unit.  We  will  explain  It  In  the  second 
section.  In  addition,  we  test  the  performance  of  various  difference  oper- 
ators with  different  post-pnocessors  on  testing  pictures. 

Fig.  2 Is  an  original  picture  without  any  disturbance  or  degradation. 
Fig.  3 Is  derived  from  Fig.  2 by  adding  an  artificial  random  noise  with  SNR  ■ 
6dB.  Fig.  If  Is  obtained  from  Fig.  2 by  replacing  the  gray  level  of  each 
point  with  the  average  gray  level  of  25  neighboring  points,  I.e., 


, 1+2  J+2 

*(I,J)  - ^ E E g(I,J) 
1-2  J-2 


Fig.  5 Is  a noisy  degraded  picture  by  adding  random  noise  with  SNR  ■ lOdB  to 
Fig.  If. 


II.  LOCAL  UNIT  TRANSFORMER 

A Local  Unit  Transformer  transforms  a picture  matrix  g(I,J)  Into 
another  picture  matrix  g'd.J)  by  the  following  formulae 

g^(>*J)  " g(i.J)/s(i,J) 


s(I,J) 


t-J+M 
K[  E 


K-I+N 

E 

KfI-N 


N,H  represent  a small  neighborhood  around  I,J  for  example  3x3  or  5x5. 
S(I,J)  actually  represents  the  energy  of  that  area  defined  by  N,M.  K Is  a 
scaling  constant. 

In  order  to  visualize  the  Idtsa  of  a local  unit  transformer,  we  can 
examine  It  In  a simple  1 dimensional  case.  Assuming  we  have  two  5 dimen- 
sional row  vectors  VI  and  V2  as: 
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A degraded  picture  from  Fig.  2 by  replacing 
the  gray  value  at  each  point  with  the  average 
of  25  adjacent  points.  I.e., 

, 1+2  J+2 


Fig.  5 A noisy  picture  from  Fig.  4 with  SNR 


VI  - [1  1 1»  10  10] 

V2  - [1  1 10  50  50] 

the  outputs  of  the  transformer  are  VI'  and  V2'  respectively.  If  we  pick  k 
as  1000,  we  will  obtain 

VI  [67.6  67.6  270  676  676] 

V2'-  [12.9  12.9  129  61»5  6^5] 

Therefore,  the  difference  between  the  maximum  and  minimum  components  of  the 
output  vectors  Is  of  about  the  same  order. 

For  further  Illustration,  we  may  consider  another  two  vectors: 

VI  - [I  1 10  1 1] 

V2  - [1  1 100  1 1]  with  K-l 

VI ' ■ r . . 1 10  10  1 . 

^vTo¥  • » AoF  » .TSTT  ^ 

V2'  ■ f ■ ^ 1 loo  1 1 . 

/I  ooo^»  /10P04  »ToooF  *Tooo¥  /lOOOA  ■* 

vr  i [0  0 1 0 0] 

V2'  i [0  0 1 0 0] 

Again,  the  difference  between  the  maximum  and  the  minimum  Is  of  the  saatc  order. 

The  examples  also  Indicate  that  the  order  of  the  applications  of  the 
transformer  and  the  difference  operator  ought  not  to  give  any  significantly 
different  results  as  far  as  edge  detection  is  concerned. 
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In  practic#  we  might  set  up  e lower  bound  (LA)  for  the  local  energy 
s(I»J).  The  transformer  Is  then  modified  as 

g'O.j)  - gO.J)/s(i,J)  If  s(lj)  > LA 

- 0 If  s(IJ)  < LA 

where  s(I,J)  Is  defined  as 

J+H  I+M  2 ,.2 

s(I,J)  - K(  Z Z gMK.A)]’'^ 

I^J-M  K-I-M 

III.  LOCAL  OPERATORS  AND  THEIR  DESICKATIONS 

In  this  section  we  define  several  local  operators  and  associate  each 
with  a proper  designation, 

(I)  Local  Unit  Transformer 

Ceflnltlori  as  In  the  previous  section. 

Designation:  T[2K<’I,  2H*’I,  K,  LA] 

(II)  Robert's  Cross  Operator 

Definition:  g^(I,j)  - |g(l-N,  j-M)  - g(l+N,  j+M)  | 

♦ |g'(»-N.  JtM)  - g(ltN.  J-M) I 

Designation:  R[2N't-|,  2H*I] 

(III)  Lap lac I an  Operator 

Definition:  g'(I J)  - |g(I-N,  j-M)  ♦ g(I-tN,  J+M) 

+ g(I-N.  J+M)  + g(l+N,  J-M)  - 4g(lJ)| 

Designation:  L[2M+I,  2H+I] 
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(iv)  E Operator 

Definition:  g'(IJ)  - [g(I,J)  - gOj)]* 

[•]*  represents  the  positive  value  of  [•]»  I.s., 

[•]*  ■ (’)  (•)<  0,  then  [•]*  ■ 0 

and  I i“I+N  k"j+M 

5n7n-- 72iff7TT2iHlT 

Designation:  E[2N+I,  2M+11 

(v)  S Operator 

Definition:  g^(I,j)  - max  (S^(l,j)  , S^ClJ) 

Sj  (I,j)  - 2Sa  - Sb  - Sc 
S2(i,j)  • 2Sa'  - Sb"  - Sc' 


If  (•)  ^ 0,  then 

gU,k) 


wher2  j+M 

Sa  - i;  g(I,i) 

1-j-M 

j+H 

Sb  - I g(l-N,i) 

il-j-M 

J+M 

Sc  - Z g(I+N,il) 

i-j-M 

j-*-N 

Sa'  - I gU,j) 

£-J-N 

j+N 

Sb'  - I g(^J-M) 

jl-j-N 

j+N 

Sc'  - z gC^-.J+M) 

t-J-N 


Designation:  S[2N+I,  2H+I] 
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I 
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(vl)  G Operator 


I 

i 


[ 


f- 

) r 


I 


I 
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Definition:  g'(l,J)  - | : 9(t,J)  - s ,(t.j)  | 

fr-I-N  j^I^I 

. J“*  J+M 

♦ I £ g(Mc)  - z g(J.k)l 
k-J-M  k-J+I  ' 


Designation:  G[2M+I,  2Mf|] 

(vll)  B Operator 

Definition:  Let  D(k)  equal  the  number  of  picture  elements  with 

gray  level  k.  P Is  a present  percentage.  ITHRESH  Is  the  smallest 
1 IMAX 

gray  level  suchthat  E w v/v  t l 

ITHRESH  ‘ ^ 

picture  element),  IMAX  Is  the  maximum  gray  lev^I. 

Then,  we  have 


g'Oj)  - I If  g(l,J)  > ITHRESH 
■ 0 otherwise. 

Dasignatlon:  B[P] 

In  addition  to  the  notations  for  operators,  we  use  PI  to  stand  for 
fig.  2,  P2,  Fig,  3;  P3,  pig,  k;  Pk,  Fig,  5,  Thus,  we  use  the  notation  B[5) 

S(3,3)  E(3,3)  PI  to  stand  for  the  picture  obtained  fay  the  following  sequence 
of  operations  on  Fig,  2,  First  we  apply  E(3,3)  on  Fig,  2 and  get  a new 
picture  Pr,  then  we  apply  S(3,3)  on  PI"  and  get  PI"",  Finally  we  apply 
3(5)  to  PI""  and  get  PI"""  , B(5]  ![3,3J  E(3,3)  P’  represents  PI"""  , 

IV,  EXPERIMENT 

A picture  a pyramid  Is  digitized  Into  a square  matrix  of  size  128  by 
128,  The  gray  level  at  each  picture  eleoenTYiSigeB  from  0 to' 2^  Thea^ pictures 
are  stored  on  magnetic  tapes.  Since  all  the  operators  except  B operator  com- 
press the  Input  picture,  we  add  some  artificial  rows  and  columns  to  the 
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picture  matrix  so  as  to  retain  the  size  of  output  picture  to  be  128x128.  For 
example,  In  the  case  of  G(5»5)  Pl»  the  top  2 rows  of  the  output  matrix  Is  made 
artificially  to  be  the  same  as  the  third  rcw;  the  lowest  two  rows  equal  the 
lowest  third  row;  the  first  two  columns  equal  the  third  column;  the  last 
two  columns  equal  the  last  third  column. 

Because  the  limitations  on  the  electrostatic  printer,  the  output  picture 
Is  rescaled  with  17  gray  levels  and  displayed  In  the  fixed-pattern  mode.  The 
experimental  results  are  shown  In  Figs.  6-35. 

V.  CONCLUSION 

Experimental  results  indicate  that  the  application  of  the  transformer 
Indeed  Improves  the  quality  of  the  output  edges  quite  a bit.  For  Instance, 

Fig.  13  and  Fig.  \k  have  almost  the  same  number  of  picture  elements,  wut 
Fig.  contains  much  more  meaningful  edges  than  Fig.  13  does.  Even  some  etges 
between  the  shaded  area  and  background (where  the  differences  In  gray  level  Is 
fairly  small  in  comparison  with  others)can  be  seen  on  Fig.  14. 

Also  It  is  found  that  the  gradient  operator  with  an  appropriate  post- 
processor can  give  a set  of  edges  which  retain  most  of  the  edge  Information 
on  the  origli^^l  picture. 

In  this  report,  we  are  dealing  mainly  with  the  behavior  of  local 
operators.  Other  questions  such  as  continuity  of  edge  points,  separation 
between  relevant  and  Irrelevant  edges,  ant  descriptions  of  edges  are  not 
considered. 
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IMt\6E  DECOHPOSITION 
J.  Burnett  and  T.  S.  Huang 

A scan  line  of  a simple  picture  Is  modeled  as  the  sum  of  two  random 
processes.  A Markov  jump  process  Is  used  to  model  abrupt  changes  In  grey 
level  due  to  boundaries.  A second  random  process  represents  variation  about 
the  average  gray  level  due  to  noise  or  texture.  The  optimum  mean  square 
filter  Is  Implementsd  to  extract  the  boundary  process.  This  estimate  of  the 
boundary  process  Is  then  differentiated  to  locate  the  jumps  or  boundary 
locations.  A comparison  Is  made  with  the  best  linear  techniques  for  perform- 
ing the  same  operations. 

Many  Image  processing  tasks  such  as  restoration,  pattern  classification 
and  data  compression  can  be  performed  with  better  results  than  otherwise  pos- 
sible If  the  Image  Is  first  segmented  or  partitioned  Into  regions  with  common 
properties.  Nahl  and  Habibi  [1]  and  Huang  and  Tang  [2j  have  suggested  this 
technique  for  Image  restoration,  Gupta  (AJ  and  Robertson  [3]  for  pattern 
classificatlai  and  Gupta  [4]  for  data  compression 

One  common  approach  to  Image  segmentation  Is  to  take  the  gradient  or 
laplacian  of  the  picture  to  find  the  boundaries.  However,  If  there  Is  much 
variation  In  gray  level  at  non  boundary  points  due  to  noise  or  texture  In  the 
Image  this  technique  does  not  work  well. 

To  overcome  this  difficulty  we  model  an  Image  l(X,Y)  as  being  the  sum  of 

two  processes:  A Markov  jump  process  (I process  whose  values  are  piece- 

wise  constant)  J(X,Y)  which  models  the  discontinuities  In  average  gray  levels 
due  to  boundaries;  and  a second  process  T(X,Y)  (where  T(X,Y)  - |(x,Y)  - j(x,Y)) 
where  T(X,Y)  represents  texture  or  noise  In  the  Image.  A sample  function  of 
a one-dimenslonal  jump  process  with  two  possible  values  Is  shown  In  Fig.  |. 
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With  this  model  the  problem  of  boundary  finding  can  be  stated  as  follows; 


given  I (X,Y)  - J(X,Y)  + T(X,Y) 

A 

what  is  the  best  estimate  J(X,Y)  of  J(X,Y)7 

A 

Once  we  have  J we  can  presumal)1y  differentiate  it  to  obtain  the  locations  of 
the  jumps  or  boundaries. 

A 

In  general  finding  J is  a very  difficult  probiem.  However,  some  one- 
dimensional results  have  been  obtained  [5]»  [6].  Let  X(t)  be  a Harkov  jump 
process  with  possible  vaiues  a^,  •••  We  assume 


(1)  P.j(h)  - P(X(t+h)  - ajlx(t)  - aj) 


where  v. . is  a constant  and 
•J 


1-Vjh  j-i 


v,jh  j»<i 


’i  ' 'ij 

j-i 

j»*l 

Let  the  observotions  be 
(2)  Y(t)  ■ X(t)  + BU(t)  where 
U(t)  is  whi te  noise. 

Let  Pj(t)  - P(X(t)  - aj|Y(S)  0 £ S<  t) 
Then  Pj  obeys  the  differentiai  equation 


dP.(t)  k 

— VjP.(t)  + -B  X [aj-X  ]Pj(t) 

+B"^[aj-X]  (3) 

where  X ■ Z aj^  P(X  I Y’(S)  0 ^ S ^ t) 

In  particular  we  consider  the  special  case  of  a random  telegraph  wave  [7]. 
This  would  be  a good  model  for  a scan  line  of  a picture  with  two  average  gray 
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levels  but  some  variation  about  each  level  due  to  texture  or  noise.  In  the 
above  formulae  we  have  k • 2,  a^  ■ I,  a2  * ~l«  v * expected  number  of  Jumps 
per  unit  tinvj  ■ I , 


X • minimum  mean  square  estimate  of  X 
- E(X(t)|Y(S)  0 < S < t)  - I«P,(t)  - I-PjU) 

From  (3)  and  (4)  we  have 

X - -2v  X -b“^X(I-X^)  + b"^(I-?)  Y(t) 


(*•) 


(5) 


A sample  of  a random  telegraph  wave  was  generated  and  Is  shown  In  Fig.  I. 

Noise  was  added  and  the  result  Is  shown  In  Fig.  2.  The  optimum  mean  square 
filter  (5)  was  Implemented  and  the  output  shown  In  Fig.  3*  The  output  of  the 
filter  was  differentiated  to  locate  the  boundaries  and  the  results  are  In 
Fig.  k.  For  comparison  the  noisy  wave  was  passed  through  a Wiener  filter 
(which  represents  the'best  linear  least  square  estimate)  and  the  output  shown 
in  Fig.  5*  Note  how  much  smoothing  of  the  signal  has  taken  place.  Tfils 
corresponds  to  the  smearing  mentioned  In  [I]  and  Is  typical  of  linear  estl** 
mators.  The  output  of  the  linear  filter  was  differentiated  to  yield  Fig.  6. 
The  boundary  spikes  of  the  optimum  filter  are  both  higher  and  narrower  than 
the  linear  filter  Indicating  that  they  would  give  more  accurate  estimates  of 
boundary  location. 

In  conclusion  we  have  shown  that  the  modeling  of  Images  with  jump 
processes  offers  the  possibility  of  developing  nonlinear  processing  techniques 
that  win  give  better  results  than  those  obtained  by  linear  methods. 
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IMAGE  SEGMENTATION  BY  UNSUPERVISED  CLUSTERING 


M.  Y.  Yoo  and  T.  S.  Huang 


I.  INTRODUCTION 

Consider  a sampled  image.  For  each  picture  element  (I.e.,  sample),  we 
measure  a set  of  n features  based  on  a neighborhood  around  that  picture 
element.  Thus,  each  picture  element  of  the  original  Image  corresponds  to  a 
point  In  the  n-dlmensionai  feature  space.  If  the  features  are  suitably 
chosen,  then  the  points  in  the  feature  space  will  form  clusters.  Unsupervissd 
clustering  techniques  can  be  used  to  partition  the  feature  space  into  severai 
regions.  Each  region,  when  mapped  back  Into  the  image  space,  gives  us  a 
(generally  disconnected)  image  segment. 

Our  study  of  the  use  of  ciustering  to  segment  images  is  partly  motivated 
by  a recently  developed  unsupervised  clustering  technique  of  Fukunaga  and 
Narendra  at  Purdue.  This  technique  has  several  advantages:  (1)  the  result  is 

independent  of  the  order  of  processing,  (ii)  the  number  of  clusters  needs  not 
be  specified,  and  (iii)  the  boundaries  of  the  clusters  lie  in  the  "valleys." 
However,  in  the  simpie  preiiminary  results  reported  beiow,  this  technique  was 
not  used. 

II.  EXPERIMENTAL  RESULTS 

In  Figure  1 we  show  an  eiectrostatic  printout  of  a digitized  image.  The 
image  consists  of  i28xi28  sampies  with  8 bit/sampie.  However,  the  electro- 
static printout  only  displays  i6  gray  levels  (using  dot-density  modulation). 
Two  features  were  measured:  the  mean  and  the  standard  deviation  of  the  3x3 
neighborhood  of  each  picture  eiement.  Figure  2 shows  the  points  in  the 
feature  space.  In  Figures  kS , some  resuits  of  ciustering  by  eyeballing  are 
shown.  In  part  (b)  of  each  of  these  Figures,  we  show  in  black  the  image 
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segment  corresponding  to  the  boxed  region  In  the  feature  space  in  part  (a). 
The  symbols  used  In  parts  (a)  of  Figures  ^-6  are  explained^  in  Figure  3- 
Ml.  FUTURE  WORK 

We  plan  to  work  on  the  segmentation  of  multispectral  images  using  the 
Fukunaga-Narendra  valley-seeking  unsupervised  clustering  technique.  We  would 
also  investigate  ways  of  incorporating  spatial  Information  (In  image  space) 
into  this  approach. 
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Standard  deviation  feature  space 


Fig.  2 Mean  - 
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NOISE  REDUCTION  IN  PHOTOGRAPHIC  IMAGES 
M.  Y.  Yoo  and  T.  S. Huang 


I.  INTRODUCTION 

Noise  reduction  using  Wiener  filters  tends  to  smear  the  edges  In  the 
image.  Here  we  report  a technique  of  noise  reduction  which  retains  edge 

sharpness. 

One  way  of  decomposing  a picture  such  that  recombination  of  the  components 
results  in  the  original  is  "synthetic  highs"  system.  In  this  system  the  image 
is  separated  into  the  "lows"  (low  frequency  part)  and  the  "edge  part"  by  the 
Gaussian  low  pass  filter  and  the  Laplacian  (or  gradient)  edge  detector.  To 
recover  the  original  the  edge  is  passed  through  the  reconstruction  filter 
(which  is  related  to  the  low  pass  filter)  before  the  recombination  with  the 
lovvs.  In  this  paper  we  use  the  synthetic  highs  system  for  reducing  the  noise 
In  the  original.  Specifically  we  used  additive  Gaussian  noise  and  generated 
noisy  picture  with  S/M  = 10  dB.  We  defined  the  slgnal-to-nolse  ratio  of  the 
image  as  the  ratio  of  the  RMS  signal  to  the  RMS  noise.  Low  pass  filter  has  a 
smoothing  effect  and  "lows"  does  not  include  much  noise.  But  the  output  of 
the  edge  r-.etector  is  quite  noisy.  If  there  is  any  way  to  remove  this  noise, 
we  c:*n  recover  a clean  picture.  The  simplest  way  is  to  threshold  the  edge 
picture.  A better  way  is  to  detect  edges  in  the  noisy  picture  and  then  use 
these  edge  locations  as  a mask  to  eliminate  noisy  points  in  the  "edge  part" 
of  the  decomposed  picture.  So  far  no  really  good  boundary  detection  algorithm 
of  noisy  picture  has  been  found.  One  fairly  satisfactory  way  proposed  by  T.»ng 
and  Huang  [l]  (see  Tang  and  Huang  in  this  report)  is  used  for  our  experiment. 
Fig.  I shows  the  schematic  diagram  of  <ur  experiment. 
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II.  LOWS 

We  used  a Gaussian  low  pass  filter  with  impulse  response 


t(x,y) 


I 


- ,Y 


2na 

The  frequency  response  Is 


2a 


r 


...  . . -2irV  (fx^  + fy^) 

(fx.fy)  - e 


If  we  define  the  bandwidth  of  the  Gaussian  filter  as  the  standard  deviation 
of  the  filter,  the  bandwidth  Is 


T 


Threshold 

Bandwidth 


Threshold  * 100 
Bandwidth  * 0.1:7^ 


Threshold 
Ban Iwidth 


Threshold  ■ 200 
Bandwidth  ■ 0.155 

1*  Reconst.  Picture  (Thresholding  used) 


5x5  Window 
Bandwidth  ' 


5x5  Window 
Bandwidth  = 0.155 


Fig.  5 Reconst.  Picture 
(Edge  detection  u 


III.  MIGMS 


The  Laplaclan  (or  gradient)  edge  detector  and  reconstruction  filter  Is 


discussed  In  Graham  [2]  and  we  just  describe  the  noise  reduction  process 


here.  The  simple  thresholding  Is  straight  forward  and  needs  no  discussion. 

Tang  and  Huang's  boundary  detection  algorithm  uses  the  combination  of 
averaging  and  the  gradient  detection.  Once  we  get  the  boundary  (from  Fig. 

2(b))  by  this  method  we  may  use  them  as  a mask  on  the  Laplaclan  picture.  All 
points  In  the  Laplaclan  picture  not  falling  In  the  boundaries  are  set  to  zero, 
tut  before  applying  this  mask,  there  Is  one  thing  to  worry  about.  Many  of 
the  pictures  have  ramp  function  type  edges  and  the  Laplaclan  edge  detector  gives 
two  spikes  (positive  and  negative)  along  the  boundary.  Therefore,  the 


boundaries  are  expanded  tc  include  both  spikes.  We  used  3x3  and  5x5 


windows  for  this  purpose:  each  boundary  point  is  expanded  to  a 3x3  ot  5x5 


block. 


The  reconstructed  picture  is  quite  improved  compared  with  the  noisy 


original.  We  also  tried  the  attenuation  of  the  high  part  before  the  recom- 


bination with  lows,  and  averaging  of  the  noisy  picture.  The  results  are 


shown  in  Figs.  2-5. 


IV.  CONCLUDING  REMAFKS 


We  used  the  Laplaclan  edge  detector  for  this  exoeriment.  But  the 


gradient  edge  detector  may  be  better  than  the  Laplaclan  for  the  noisy  pic- 


ture and  we  will  continue  this  work  with  the  gradient. 


REFERENCES 


[1]  G.  Y.  Tang  and  T.  S.  Huang,  "Edge  Extraction  Techniques,"  this  Report, 


[2]  D.  N.  Graham,  "Image  Transmission  by  Two-Dimensional  Contour  Coding," 
Proc.  IEEE,  Vol.  55,  No.  3,  March  1967. 


P 


92 


SEGMENTING  COLOR  PICTURES  WITH  BLOB 


Paul  A.  Wintz 

The  standard  color  pictures  shown  In  Fig.  1 was  first  transformed  from 
the  red,  green  and  blue  coordinate  system  through  Y,  i,  and  0 coordinate 
system.  The  BLOB  algorithm  was  then  run  on  the  Y,  1,0.  pictures.  All  of  the 
picture  elements  in  each  BLOB  in  the  Y,  I,  Q dc>main  were  replaced  with  an 
average  value.  The  Inverse  transform  from  the  Y,  i , Q domain  back  to  the  red, 
green  and  blue  domain  was  then  made.  The  resulting  pictures  for  several 
different  settings  of  the  thresholds  are  represented  in  Figs.  2,  3»  snd  4. 

The  reconstructed  pictures  in  Figs.  2,  3*  and  k show  only  a smiili  amount  of 
distortion  from  the  original  indicating  the  possibility  of  using  BLOB  as  a 
coding  algorithm  since  significant  compression  ratios  can  be  obtained  by 
coding  the  BLOB's  relative  to  coding  the  original  picture  elements 

Reference 

1.  P.  A.  Wintz  and  J.  N.  Gupta,  "A  Boundary  Finding  Algorithm  and  its 
Applications,"  !-EE  Trans,  on  Circuits  and  Systems,  March  1975. 


STATISTICAL  MEASURES  FOR  THE  CLASSIFICATION  OF  TEXTURES 
0.  R.  Mitchell  and  C.  R.  Myers 

The  need  to  classify  Image  regions  has  led  us  to  search  for  optimum 
statistical  measures  (or  features)  for  the  classification  of  Image  textures. 
First,  we  have  sought  to  Improve  existing  measures  and  second,  we  have 
derived  promising  new  measures. 

The  statistical  measures  of  texture  based  on  spatial  dependence  probabil- 
ities (or  the  neighboring  grey  tone  co-occurence  properties)  has  been  one  of 
the  most  promising  [1].  However,  these  measures  are  not  Invariant  to  resolu- 
tion and  most  of  them  are  not  Invariant  to  Illumination  changes. 

BRIGHTNESS  INVARIANCE 

The  effect  of  a change  In  Illumination  level  (Including  shade  vs  sun  and 
film  development  differences)  Is  primarily  a multiplicative  effect.  To  make 
our  texture  measures  Invariant  to  this  effect  we  have  taken  the  logarithm  of 
all  Intensities  and  subtracted  the  local  mean.  This  Is  a form  of  homomorphic 
filtering  and  Is  essentially  equivalent  to  the  first  processing  stages  In  the 
human  visual  system  (a  logarithmic  transformation  followed  by  a band-pass 
filter)  [2]. 

For  Initial  tests  we  have  used  texture  samples  from  a book  by  Brodatz  [3l, 
Figure  1 shows  typical  textures.  We  used  ten  different  textures,  5 samples 
of  each  for  training  and  2 samples  of  each  for  testing.  Each  sample  was 
bkxbk  points,  with  6k  grey  levels.  The  features  used  for  classification 
were  those  suggested  by  Haralick  [2].  The  classification  procedure  was  a 
simple  normalized  Euclidian  distance  measure  with  all  features  weighted 
equally.  The  results  are  shown  In  Fig.  2.  The  logarithmic  measure  gave 
results  comparable  to  the  standard  method  with  these  textures.  However,  the 
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RESOLUTION  INVARIANCE 


In  order  to  make  the  spatial  dependence  matrix  measures  of  texture  Invar- 
iant to  resolution  (Includes  camera  or  scanner  distances  and  sampling  rate)  we 
have  used  the  two-dimensional  regional  autocorrelation  to  regulate  the  dis- 
tances used  In  calculating  the  spatial  dependence  matrix.  Typical  autocor- 
relations (actually  autocovariances)  are  shown  In  Figure  3»  We  have  used  the 
contour  at  which  the  autocorrelation  has  fallen  to  a fixed  percentage  of  its 
maximum  value  as  the  primitive  texture  shape.  This  contour  Is  then  used  for 
selecting  horizontal  and  vertical  distances  to  be  used  In  calculating  the 
spatial  dependence  matrix.  Results  using  the  same  texture  training  and  test 
samples  as  before  are  shown  In  Figure  k.  Vfhen  a contour  of  20^  maximum  value 
Is  used  the  results  are  very  encouraging.  This  measure  Is  Invariant  to  changes 
In  resolution  and  sampling  rate  (as  long  as  the  Image  Is  not  undersampled). 
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Fig.  k Variabie  Distance  Spatial  Dependence 
Classification  Results 


A NEW  STATISTICAL  MEASURE 

A new  measure  for  texture  classification  has  been  developed  which  is 
computationally  much  simpler  than  the  above  methods.  It  is  based  on  the  human 
visual  system  Intuition  that  the  Important  texture  information  Is  contained  In 
the  relative  frequency  of  local  extremes  of  various  sizes  In  Intensity.  The 
principal  measure  In  this  process  is  the  determination  of  the  number  of  max- 
ima and  minima,  along  a one-dimensional  scan  direction.  The  extremes  are 
not  included  unless  tliey  meet  a threshold  condition  as  shown  in  Figure  5. 

The  maximum  is  called  a maximum  only  If  the  Intensity  falls  the  threshold 
amount  below  the  maximum  before  a higher  valued  Intensity  Is  encountered. 

This  hysteresis  effect  Ignores  small  ripples  while  including  the  larger  ones. 
By  repeating  the  process  for  several  threshold  settings,  a vector  of  numbers 
Is  obtained  which  characterizes  the  texture. 

We  have  found  that  the  steps  necessary  to  make  this  measure  Invariant 
to  liiumination  and  resolution  are  very  easily  Implementable  and  improve  the 
accuracy  of  the  results.  These  steps  are  (I ) take  the  logarithm  of  the 
Intensities  first;  (2)  use  the  ratios  of  the  number  of  extrema  at  each 
threshold  to  the  next  Instead  of  numbers  of  extrema  themselves.  The  threshold 
levels  used  for  these  measures  are  shown  In  Figure  6.  The  results  are  shown 
In  Figure  7.  We  used  7 features  In  both  horizontal  and  vertical  directions 
for  a total  of  features.  The  results  are  very  encouraging  indeed. 
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iVlAX.  - MIN.  SELECTION 
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■ MIN.  DETECTED 


Fig.  5 MAX-MIN  Selection 
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CONCLUSION 


We  have  derived  feature  classification  techniques  which  are  Invariant  to 
Illumination  and  resolution  changes.  A new  technique  Is  described  which 
allows  simple  computation  of  texture  features.  Although  our  sample  sizes 
are  small,  we  feel  the  results  do  show  that  the  techniques  described  can  be 
optimized  to  provide  large  Improvements  In  texture  classification  techniques. 
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AN  OPTIMAL  FEATURE  SUBSET  SELECTION  ALGORITHM 
K,  Fukunaga  and  P.  M.  Narendra 


Patterns  to  be  classified  by  a recognition  system  are  generally  chara- 
cterized by  a set  of  measurements  or  features.  Often,  the  dimensionality 
of  this  feature  space  Is  too  large  for  efficient  and  reliable  application  of 
existing  classification  techniques.  The  feature  extraction  problem  is  then, 
to  reduce  the  dimensionality  of  the  feature  space,  without  significantly 
affecting  the  discriminatory  capacity  of  the  feature  set.  Classification 
with  the  reduced  set  of  features  is  then  facilitated. 

One  approach  to  feature  extraction  is  to  select  a smaller  subset  of  m 
features  of  the  set  of  n original  features  (m  < n).  There  are  ml  (n-m) ! 


such  subsets  of  size  m.  The  best  subset  of  size  m is  then  the  subset  for 

which  the  value  of  a criterion  is  extremized  over  all  subsets  of  size  m. 

One  way  to  obtain  the  best  subset  is  of  course,  to  evaluate  the  criterion  for 

ali  the  (")  subsets  and  choose  the  one  that  yields  the  extreme  value.  But 
m 

exhaustive  search  is  computationally  unfeasible  even  for  modest  values  of 
n and  m,  as  the  amount  of  search  involved  becomes  prohibitive.  For  example, 
with  n*2l*  and  m-12,  the  number  of  subsets  to  be  considered  is  2,70^,156. 
Several  approaches  have  been  suggested  to  select  a good  subset  without 
exhaustive  enumeration.  Forward  selection  [i],  and  Dynamic  Programming 
techniques  [2,3]  are  among  these.  But  none  of  these  techniques  guarantee 
that  the  selected  subset  is  indeed  the  best  among  all  subsets  of  size  m. 

We  have  formulated  the  feature  subset  selection  problen  as  a combina- 
tional optimization  problem.  A branch  and  bound  search  algtrithm  [A,5l  has 
been  developed  to  find  the  best  subset  without  exhaustive  search.  The 
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algorithm  is  efficient  and  is  guaranteed  to  be  optimal.  The  algorithm  in- 
volves recursive  enumeration  of  the  subsets.  Unfeasible  subsets  are  dis- 
carded without  being  explicitly  enumerated.  The  algorithm  is  stated  for  any 
general  criterion;  recursive  equations  are  deveioped  to  aid  the  implementa- 
tion of  the  aigorithm  for  a class  of  quadratic  criteria  such  as  the  Fisher 
criterion,  Divergence  and  the  Bhattachar^’ya  Distance. 

To  test  the  effectiveness  of  the  fe2ture  selection  algorithm,  it  was 
applied  to  mul tispectral  data  from  airborne  scanners.  The  data  originated 
from  the  Laboratory  for  Applications  of  Remote  Sensing  (LARS)  and  consisted 
of  classified  data  points  of  fieids  pianted  with  corn,  soybean,  rye,  wheat, 
etc.  The  problem  was  to  select  out  of  the  12  data  channels  encompassing  the 
entire  spectrum  (0.4  to  1.0  micrometer)  the  subset  of  4 channels  that  yields 
the  optimum  values  of  a criterion.  The  Fisher  criterion  was  chosen,  Exhaus- 
tive  search  involves  the  computation  of  ail  the  ( ^^  ) ■ 495  subsets  of  channels 
to  obtain  the  best  subset.  The  present  aigorithm  was  applied  to  obtain  the 
best  subset  to  di scrim! nate  soybean  from  corn.  The  computational  effort  was 
reduced  by  a factor  of  more  than  50  over  exhaustive  enumeration,  rt  the 
same  time  guaranteeing  that  the  chosen  best  subset  w^s  indeed  the  best  among 
ail  the  495  subsets  of  size  4. 

Further,  the  aigorithm  was  applied  to  choose  the  best  subset  of  12  out  of 
24  channels.  There  are  2,704,150  candidate  solutions.  The  algorithm 
obtained  the  best  subset  with  the  computational  effort  equivalent  to  com- 
puting the  criterion  for  6000  subsets.  The  savings  are  indeed  substantial. 

This  work  is  being  documented  for  publication.  Further  work  is  underway 
to  extend  the  approach  to  other  problems  In  Pattern  Recognition. 
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EXTENSION  OF  THE  SUBSPACE  METHOD  IN  PATTERN  RECOGNITION  TO  PICTURES 

G.  V.  Sherman  and  K.  Fukunaga 
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I.  INTRODUCTION 

Most  decision-theoretic  pattern  recognition  schemes  are  of  the  zonal  va- 
riety. That  Is,  observations  are  made  In  an  n-dimenslonal  vector  space  which 
has  been  partitioned  Into  M zones,  corresponding  to  the  M classes,  by  (n-l)- 
dlmenslonal  surfaces.  Fig,  1 Illustrates  this  concept  for  n«2  and  M“3» 

Ideally  the  partitions  In  Fig.  1 are  determined  by  Bayes  Law.  However, 
usually  the  n-dlmensIonal  class  dependent  densities  are  unavailable  and  too 
cutnbersome  to  compute.  Hence,  suboptimum  methods  are  used  to  determine  the 
partitions. 

Watanabe  [1]  has  proposed  a different  approach.  Instead  of  partitioning 
the  observation  space  Into  zones,  perhaps  one  should  find  feature  subspaces 
which  correspond  In  a one-to-one  manner  with  the  classes.  Fig.  2 Illustrates 
this  concept  for  n-3  and  M*2. 

Wantanabe's  original  approach  can  be  described  as  follows: 

1.  Construct  a representation  subspace  for  each  class  from  the  eigenvectors 
of  the  class  autocorrelation  matrices;  I.e.,  the  discrete  Karhunen-Loeve, 
(K-L)  expansion.  Retain  only  these  eigenvectors,  corresponding  to  the 
largest  eigenvalues,  which  contain  t%  of  the  sample  energy  on  the  average. 
Typically  t ~ 95  and  Is  the  same  for  each  class. 

2.  Construct  a feature  subspace  for  each  class  by  subtracting  from  the  cor- 
responding representation  space  Its  Intersections  with  the  representation 
spaces  of  other  classes.  The  feature  space  are  then  transformed  so  that 
they  are  all  mutually  orthogonal. 

3.  That  portion  of  the  observation  space  which  Is  not  Included  In  any 
feature  space  makes  up  the  reject  space. 
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4.  Unknown  samples  are  projected  onto  the  various  class  feature  spaces  and 
classl f Icotlon  Is  accomplished  by  a Bayes  rule  which  Is  based  on  member- 
ship In  a subspace  rather  than  membership  In  a zone. 

Therrlen  [2]  made  the  subspace  method  more  practical  by  deriving  effi- 
cient method  for  computing  the  intersection  and  union  of  subspace. 

II.  EXTENSION  TO  PICTURES 

The  representation  spaces  for  different  classes  of  pictures  cannot  be 
obtained  by  the  Karhunen-Loeve  approach.  The  dlnenslon..!  1 t'.es  Involved  are 
too  high.  Therefore,  the  pictures  are  represented  by  an  outer-product 
expansion  of  the  form 

n n _ 

A - Z E Y, , <I>:  (0 

1-1  J-1  J 


where  A Is  an  nxn  raster-scanned  picture,  y.j  Is  a scalar;  and  are 

n-d  1 mens  Iona  1 vectors.  The  { (Jij  } are  { ifij  } are  complete  orthonormal  bases  which 

span  the  column  and  row  spaces  of  A respectively.  That  Is,  any  column  of  A 

can  be  expressed  as  a linear  combination  of  (J)^'s  and  similarly  for  rows  of  A 

and  The  { } and  { t{)j  } bases  can  be  obtained  by  many  methods,  fixed 

bases  such  as  Fourier  and  Hadamard,  or  data  dependent  bases  such  as  discrete 

2 

K-L.  The  latter  was  choosen  for  our  experiments.  Note  than  an  n -dimensional 

K-L  expansion  has  been  replaced  by  two  n-dlmenslonal  K-L  expansions. 

Furthermore,  we  have  shown  that  the  Intersection  and  union  of  outer- 

product  subspaces  can  be  calculated  by  finding  the  Intersection  and  union  o^ 

the  corresponding  column  and  row  spaces  separately.  Again  the  result  Is  that 
2 

an  n -dimensional  eigenvector  calculation  has  been  replaced  by  two  n-dlmen- 

sional  eigenvector  calculations.  This  result  does  not  depend  on  the  choice 

of  { <{»,  } and  { lb.  }. 

• J 
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Featurt:  extraction  Is  completed  by  transforming  the  feature  opaces  so 


that  they  are  all  mutually  orthogonal.  This  operation  originally  would 
2 

require  an  n -dimensional  matrix  Inversion.  But  we  have  shown  It  to  be 
equivalent  to  two  n-dimenslonal  matrix  Inversions  operating  again  on  the  row 
and  column  spaces  separately. 

The  subspace  classifier  proposed  >'y  Watanabe  [I]  was  based  on  heuristic 
arguments.  It  requires  the  calculation  of  the  following  quantities: 

1.  p(A|u)|^),  the  probability  of  a class-k  sample  "appearing"  In  subspace  Jl. 

2.  p(tjJj^),  the  a priori  probability  of  class  k. 

3.  q((U|^|£),  the  probability  of  a sample  being  In  class  given  that  It 
"appears"  In  subspace  £. 

k.  W(i|)^),  the  probability  of  random  sample  appearing  In  subspace  SL, 

5.  <1  |25) » the  probability  of  X_  being  from  class  coj^. 

A serious  question  arises  concerning  this  classifier.  Can  we  define  the 
event  ")^  appears  In  subspace  £"?  In  most  cases,  ur  known  sample  )^wlll  not  be 
completely  contained  In  any  feature  space.  However,  we  Intuitively  feel  that 
appears  In  subspace  to  a greater  extent  than  it  appears  In  subspace 
If  the  projection  of  £ on  has  a larger  norm  than  that  of  the  projection 
of  £ on  S^.  Thus  ")^  appears  In  subspace  il"  Is  a fuzzy  event  [3,^]»  Zadeh  shows 
how  to  define  probability  of  a fuzzy  event  and  also  defines  conditional 
probability  and  Bayes  rule  for  fuzzy  events.  This  formalism  results  In  a 
different  definition  for  the  empirical  calculation  of  p(2|W|^).  All  other 
quantities  remain  the  same. 

III.  EXPERIMENTAL  RESULTS 

The  pictures  used  for  verification  of  the  above  theory  were  the 
numerals  from  the  Munson  alphanumeric  data  available  through  Stanford 
Research  Institute.  Since  we  are  not  Interested  In  Implementing  a 
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character  reader,  only  the  numerals  were  used.  Each  numeral  Is  a 24x24  binary 
matrix.  147  samples  of  each  numeral  are  available.  98  samples  were  used  for 
design  and  4S  for  testing. 

Wantanabe's  method  for  finding  feature  subspaces  was  found  to  be 
computationally  unfeasible,  and  Therrlen's  method  had  a theoretical  flaw. 
Therrlen  proposed  that  the  1th  feature  space  S'  should  be  obtained  from  the 


representation  space  by 


M 


S'  - S,  ( n s ) 
' k-1 

kj<l 


(2) 


where  Is  the  complement  of  S^.  This  approach  resulted  In  empty  featuro 
spaces  because  subspaces  do  not  obey  the  distributed  law  of  logic  given  by 


That  Is  Sj  n - 0 does  not  Imply  Sj  n ^1^  Therefore,  a layered 

machine  was  Implemented  In  which  feature  subspaces  are  formed  by 


(3) 


Si  ns^ 


(4) 


Tha*.  Is,  the  problem  was  reformulated  Into  the  form  shown  In  Figure  3. 

For  our  10-class  problem,  the  layered  approach  reduced  the  number  of 
Intersection  calculations  from  90  to  9.  Furthermore,  (4)  resulted  In  non- 
empty feature  spaces.  The  classifier  Itself  Is  yet  to  be  Implemented. 

IV.  FUTURE  WORK 

As  mentioned  above,  the  features  must  be  evaluated  by  a subspace 
classifier.  Moreover,  It  Is  believed  that  these  experiments  can  reveal  the 
relationship  between  observation  space  dimensionality,  number  of  classes. 
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V ,,  p 
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and  number  of  training  samples  necessary  for  valUi  classifier  design  and 

valid  feature  extraction. 
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FOURIER  DESCRIPTORS 
Paul  Wlntz  and  Lyle  Stanfield 

Shape  is  an  important  feature  for  pattern  recognition  algorithms.  One 
promising  method  of  shape  description  Is  to  code  a contour  (outline  of  a 
shape)  as  a sequence  of  complex  numbers  and  to  work  with  the  Fourier  trans- 
form of  this  sequence. 

The  descriptors  used  In  this  project  are  identical  to  those  described 
by  Granlund  [l].  A contour  is  plotted  In  the  complex  plane.  As  the  contour 
is  traced  a complex  function  of  arc  length  is  obtained  from  the  coordinates 
of  the  points  on  the  contour.  The  function  Is  periodic  since  it  Is  pos- 
sible to  trace  around  the  contour  more  than  once.  The  function  Is  normalized 
so  that  its  period  (or  the  length  of  the  contour)  is  2ir  . The  Fourier  trans- 

form is  taken  and  is  known  as  the  Fourier  descriptor  of  the  shape. 

The  shapes  studied  in  this  project  were  typed  on  cards  and  read  into  a 
50x50  array  in  a FORTRAN  program.  A routine  was  written  to  find  an  initial 
point  and  to  trace  the  contour.  The  coordinates  of  the  points  on  the  contour 
were  stored  in  a complex  array  whose  length  was  adjusted  to  be  a power  of  two. 
The  adjustment  was  required  by  the  Fourier  transform  routine  and  was  accom- 
plished by  assuming  that  points  on  the  contour  were  connected  by  straight  line 
and  interpolating  between  points.  The  interpolation  caused  some  rounding  of 
the  corners  of  objects. 

The  magnitude  and  phase  of  the  Foruler  coefficients  were  plotted  on  the 
Gould  plotter.  The  magnitude  plots  are  log-log.  The  a^  coefficients  are 
shown  as  *'s  on  the  y-axes. 

A routine  was  also  written  to  set  all  but  certain  coefficients  to  zero. 

The  Inverse  transforms  were  then  taken  to  Identify  the  shapes  causing  certain 
'envelopes'  In  the  coefficient  plots. 
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Tile  shapes  checked  were  square,  cross,  circle,  and  triangle.  The 

circle  and  triangle  were  distorted  because  of  the  coarse  50x50  array  size. 

The  initial  point  of  aii  contours  was  the  left-most  point  at  the  bottom  of 

the  contour.  Contours  were  traced  in  a counter-clockwise  direction. 

Fig.  1 shows  the  square  as  it  was  before  taking  the  transform.  Fig.  2 

shows  the  magnitude  and  phase  of  the  Fourier  descriptor.  On  the  magnitude 

plot  note  that  there  are  four  distinct  envelopes.  Every  fourth  point  seems 

to  be  related.  Also  note  that  the  coefficients  a where  n ■ 1 + 4p  are  much 

n 

larger  than  the  other  coefficients,  as  required  by  Graniund's  symmetry  con- 
dition. 

Figs.  3”6  show  reconstructions  of  parts  of  the  square.  Fig.  3 is  the 


inverse  transform  of  the  coefficients  a^  a^,  a^ a^+^*k.  Aii  other  coef- 


ficients were  set  to  zera.  Due  to  the  square's  symmetry  these  coefficients 
should  contain  aii  of  the  information  about  the  square.  Notice  that  aii  of  the 
corners  are  rounded  identically  so  the  contour  has  a greater  amount  of  sym- 
metry than  the  original  square. 

Figs.  ^*-6  were  generated  by  the  sets  of  coefficients  of  the  forms 
^2+l*k'  ^3+i*k’  1 L/^»).  Note  that  the  contours  have  rotation- 

al symmetry,  but  they  don't  meet  Graniund's  necessary  condition  for  rotation- 
ally  symmetric  objects.  This  disc-epancy  appears  frequently  in  the  remaining 
plots  and  wi i i be  explained  later. 

Figs.  7"i2  are  for  a cross-shaped  contour.  They  correspond  exactly  to 
the  plots  for  the  squares  and  have  many  of  the  same  properties. 
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There  are  two  reasons  why  some  of  the  generated  shapes  had  rotational 
symmetry  without  meeting  Graniund's  condition  for  symmetry.  The  first  reason 
is  that  Granlund  assumed  that  the  contour  was  traced  only  once.  He  also 
assumed  that  the  lines  of  the  con^’ours  did  not  cross  over  themselves  as  in  the  star. 
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I These  conditions  are  met  when  descriptors  are  calculated  from  a contour,  but 

they  may  not  be  met  when  a'contour  Is  generated  from  an  artificial  set  of 
descriptors. 

A simple  example  of  a contour  being  multiply  traced  would  be  a circle 
generated  from  a Fourier  component  other  than  aj.  The  parametric  represen- 
tation of  a circle  consists  of  sine  and  cosine  components  with  equal  ampli- 
tudes. The  component  a^  will  generate  a circle  which  will  trace  once  around  Its 
contour  in  a period  of  2tt.  The  coefficient  a^  has  twice  as  high  a frequency, 
so  it  will  trace  around  a circle  twice  in  the  period  2ir.  The  resulting  cir- 
cles look  alike  even  though  they  have  different  coefficients  and.  In  a certain 
sense,  different  contours. 

The  properties  of  contours  generated  by  descriptors  whose  only  non-zero 

coefficients  are  of  the  form  a.^  will  now  be  examined.  The  symbols  stand 

i+pm 

I for: 

i jnitial  component,  lowest  frequency  component 
m multiples,  every  mth  component  may  be  non-zero 
p integer  variable 

The  following  symbols  will  be  useful: 

s degree  of  rotation  symmetry  of  the  contour 
z the  number  of  times  the  contour  is  traced  in 
a period  of  27T 

I c the  number  of  nonoverlapping  cycles  made 

! around  the  center  of  the  contour 

This  derivation  is  based  on  Granlund's  idea  of  rotating  the  object  and 

moving  the  starting  point  so  that  the  effects  of  the  two  changes  cancel  out. 

Let  the  initial  point  on  a rotational ly  symmetric  contour  be  moved  to  the 

next  corresponding  point  along  the  contour.  2ir/z  is  the  distance  to  move  to 

trace  once  completely  around  the  contour,  so  a move  of  2tt/zs  will  result  in  a 

move  to  the  next  point  which  Is  similar  to  the  original  initial  point.  Note 

that  this  move  is  to  the  next  similar  point  on  the  contour.  This  may  not  be 
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the  point  at  an  angle  2tt/s  radians  back  around  the  contour.  For  example,  with 


Rotating  the  contour  by  -2ttc/s  radians  will  move  the  Initial  point 
back  to  Its  original  location  In  the  complex  plane  and  will  result  In  a shape 
with  the  same  Fourier  descriptors  as  the  original. 


For  this  equation  to  hold  either  a^  - 0,  or  (n/zs)  - (c/s)  ■ for  Integer 
D > 0.  Solvinq  for  n gives  n - zc  + pzs.  This  result  Is  of  the  form  n ■ I; 


which  Is  the  type  of  Fourier  descriptor  we  are  Interested  In. 

The  definition  of  c requires  that  It  not  Include  the  result  of  multiple 
traces  (traces  which  write  directly  over  previous  traces).  If  c and  s were 
not  relatively  prime  then  multiple  tracing  would  result  before  the  first  c 
cycles  were  complete.  Therefore  c and  s must  be  relatively  prime.  Therefore, 
the  greatest  common  divisor  of  i and  m Is  z. 

Relations  which  must  hold  for  rotational ly  symmetric  contours 


number  of  tracings 
degree  of  symmetry 
number  of  distinct  cycles 
total  number  of  cycles 


These  results  show  how  to  calculate  the  degree  of  rotational  symmetry, 
number  of  cycles  about  the  center  of  the  contour,  and  the  number  of  multiple 


traces  around  the  contour,  given  that  the  contour  is  symmetric.  Unfortunctely 
there  are  no  known  conditions  which  are  sufficient  to  require  a generated  con- 
tour to  be  symmetric. 

These  results  reduce  to  Granlund's  result  when  2 anti  c are  both  1.  The 
equations  above  explain  the  synmetry  of  all  of  the  test  cases  run  for  this 
project.  The  contours  generated  by  the  circle's  descriptors  are  especially 
interesting  because  of  the  range  of  their  degrees  of  symmetry. 
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Fig.  i Contour  of  the  input  data 
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Fig.  3 The  contour  generated  by  and  coefficients  of  the  form  *^44^1^' 
where  K Is  an  Integer.  ' The  other  coefficients  t#ere  set  equal 
to  zero. 
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Fig.  5 The  contour  generated  by  and  coefficients  of  the  form 

where  K Is  an  Integer.  The  other  coefficients  were  set  equal 
to  zero. 
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FREaUENC  : IN  R«JIflH5/5EC 
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Fig.  9 The  contour  generated  by  and  coefficients  of  the  form 

where  K i s an  Integer.  The  other  coefficients  were  set  equal 
to  zero. 
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Ffg.  10  The  contour  generated  by  and  coefficients  of  the  form  A2+4K 
where  K Is  an  Integer.  The  other  coefficients  were  set  equal 


to  zero. 
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Fig.  II  The  contour  generated  by  and  coefficients  of  the  form 

where  K I s an  Integer.  The  other  coefficients  were  set  cquaT 
to  zero. 
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Fig.  12  The  contour  generated  by  a and  coefficients  of  the  form  a,  . 

O 

where  K Is  an  Integer.  The  other  coefficients  were  set  equal 
to  zero. 


138 


SELECTION  OF  SPECTRAL  BANDS  FOR  USE  IN  MULT  I SPECTRAL  SCENE  ANALYSIS 

E.  WI swell  and  G.  Cooper 


3 


I.  INTRODUCTION 

Selection  of  spectral  bands  In  multlspectral  scanners  has  to  the  present 
time  oeen  based  largely  on  available  energy  considerations.  Hence,  selec- 
tion of  spectral  bands  with  consideration  to  numerical  scene  analysis  has  to 
date  been  largely  empirical.  Therefore,  this  investigation  Is  directed  to- 
ward selection  of  bands  In  multlspectral  scanners  quantitatively  with  regard 
to  the  use  of  the  data  for  numerical  processing.  To  optimally  select 
spectral  bands,  a criteria  of  optimality  must  accordingly  be  chosen  and 
evaluated.  Thus,  this  project  Is  concerned  with  optimal  multlspectral  band 
selection  with  regards  to  numerical  scene  analysis  and  a suitable  optimality 
criteria. 

II.  AREA  OF  INVESTIGATION 

Presently,  this  investigation  Is  attempting  to  properly  apply  Infor- 
mation theoretic  concepts  to  the  problem  of  multlspectral  band  selection. 

A first  step  Is  to  consider  the  spectral  response  as  a random  process  In 
wavelength.  Hence,  represent  the  spectral  response  as: 

y(X)  - s(X)  + n(X) 

Here  s(X)  Is  the  true  spectral  response.  n(X)  Is  a term  that  Includes  noise 
(i.e.  thermal  noise)  as  well  as  any  Interference.  An  example  of  Interference 
might  be  bare  soil  between  rows  of  crop  vegetation  early  In  the  growing  sea- 
son. It  Is  Important  to  carefully  consider  the  consequences  of  declaring 
part  of  the  received  spectral  response  as  noise  when  It  may  actually  contain 
Information  useful  for  numerical  scene  analysis.  Perhaps  It  would  be  useful 
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to  decompose  n(X)  Into  a noise  term  and  an  Interference  term.  This  thought 
remains  to  be  Investigated, 

A spectral  interval  (X^.X^)  where  X Is  wavelength  Is  considered.  The 
Interval  Is  to  be  subdivided  Into  several  spectral  bands.  It  Is  proposed  to 
compute  the  mutual  Information  between  spectral  bands  and  to  choose  a subset 
of  spectral  bands  which  maximizes  the  mutual  Information  between  the  remaining 

bands. 

In  order  to  compute  the  mutual  Information  between  spectral  bands,  the 
autocorrelation  functions  of  the  signal  R^CX^.X^)  and  R„(X,,X^)  (corresponding 
to  the  n(X)  term)  are  needed.  Therefore,  an  effort  Is  currently  under  way  to 
estimate  the  autocorrelation  functions  from  representative  mul tispectral  data 
obtained  at  Purdue's  Laboratory  for  Applications  of  Remote  Sensing  (LARS). 

A problem  related  to  the  computation  of  mutual  Information  Is  consid- 
eration of  the  effects  of  width  of  the  spectral  bands.  Another  related  prob- 
lem concerns  the  relationship  of  the  spectral  bands  to  temporal  and  spatial 
Information.  The  effects  of  a set  of  mul tispectral  bands  and  their  width 
will  ultimately  be  related  to  spatial  resolution  for  example.  And,  of  course, 
there  Is  no  reason  to  suppose  that  spectral  Information  will  be  constant 

temporally. 

III.  FUTURE  INVESTIGATION 

There  are  two  main  topics  to  be  pursued.  First,  the  Information 
theoretic  approach  will  be  continued.  Secondly,  It  may  be  useful  to  pursue 
a decision  theoretic  approach.  The  motivation  for  this  that  one  may  be  able 
to  choose  a cost  function  that  reflects  the  application  for  which  the  data  Is 
to  be  used  (I.e.,  classification).  The  extent  to  which  any  decision 
theoretic  approach  differs  from  previous  techniques  (I.e.,  pattern  recogni- 
tion approaches)  must  be  determined. 
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IV.  COMMENTS 


r- 

‘ Due  to  the  fact  that  this  Investigation  Is  In  Its  Initial  stages, 

w application  oriented  results  have  yet  to  be  achieved.  However,  a theoretical 

approach  that  Is  perceived  as  viable  has  been  formulated.  The  experimental 
demonstration  of  this  approach  Is  to  be  pursued. 


ITERATIVE  IMAGE  RESTORATION,  III 
T.  S.  Huang,  S.  Berger,  and  M.  Kaveh 


. INTRODUCTION 


The  restoration  of  Images  degraded  by  linear  space''var lent  degradations 
can  be  accomplished  through  the  use  of  the  projection  algorithm.  This 
algorithm  has  been  applied  to  both  one  and  two-dimensional  signals  In  order 
to  examine  Its  capabilities.  The  current  results  Indicate  that  the  method 
may  be  quite  valuable  In  the  restoration  or  Images. 

1 I . INVESTIGATION  OF  THE  ALGORITHM 

The  projection  algorithm  is  basically  an  Iterative  method  for  the 
solution  of  a system  of  equations 


g,  - f,  fj  + ...  *8,^ 


■ *21  ^1  * ••• 


♦ a-  f 
2n  n 


3m'  *ml  ^ * ••• 


+ 5 f 
mn  n 


where  n may  not  be  equal  to  m.  The  solution  Is  obtained  by  successive  pro- 
jections onto  planes  In  hyperspace.  For  the  case  ' "'ere  a unique  solution 
does  exist,  the  algorithm  will  converge  to  the  point  of  Intersection  of  the 
hyperplanes.  If  the  planes  do  not  Intersect  at  a single  point,  the  algorithm 
will  converge  to  a point  which  may  be  a useful  approximation  In  a restoration 


sense. 


In  order  to  apply  the  algorithm,  one  must  consider  the  degradation  pro- 
cess In  terms  of  discrete  components.  A one-dimenslonal  signal  Is  divided 


|l*2 


into  eveniy  spaced  samples,  ard  a two-dimensional  Image  Is  treated  as  an 
array  of  picture  elements.  The  degradation  process  consists  of 
a known  linear  degradation  and  additive  Gaussian  noise. 

The  performance  of  the  algorithm  was  Investigated  for  several  cases.  A 
one-d Imens Iona  1 signal  was  degraded  by  a low-pass  filter.  The  original 
signal  f(x)  consists  of  two  unit  pulses  of  five  points  each,  separated  by  two 
zero  points  (Fig.  1).  The  degradation  is  effected  in  the  Fourier  transform 
domain.  The  discrete  Fourier  transform  (OFT)  of  f(x)  Is  obtained  by  a fast 
Fourier  transform  routine.  A triangular  low-pass  filtering  operation  (with 
a cutoff  frequency  equal  to  .22  times  the  maximum  frequency)  yields  the 
degraded  signal  gj (x)  (Fig.  2).  The  result  after  the  addition  of  Gaussian 
noise  Is  denoted  by  g(x)  (Fig.  3).  The  Initial  guess  for  the  algorithm  Is 
taken  to  be  g(x).  The  algorithm  was  tested  for  different  restrictions  on  the 
successive  Iterations  f*(x). 

The  projection  algorithm  was  compared  with  the  least-squares  Inverse 
filter,  whose  frequency  response  Is  given  by 

H*(u)  • Sr(u) 

Q(u)  3 

|H(u)r  • S^(u)  + S^(u) 

where  H(u)  Is  the  degradation  filter,  S^(u)  and  S^(u)  are  the  spectral 
densities  of  the  noise  and  signal,  respectively. 

A two-dimensional  signal  was  also  used  to  test  the  projection  algorithm. 
An  image  was  degraded  by  an  averaging  operation.  Each  picture  element  was 
represented  by  an  Integer  from  0 to  255»  which  Is  8 bit  quantization.  The 
128x128  image  arrays  were  stored  on  magnetic  tape  In  a compacted  format.  The 
Gould  electrostatic  plotter  was  used  to  obtain  half-tone  reproduction  of  the 
images. 
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F H . RESULTS 


The  restoratFon  was  FmpFemented  on  a CDC  6500  computer.  For  the  one- 

dFmensFonaF  case,  the  addFtFve  GaussFan  noFse  had  a standard  devFatFon  of 

^ ■ 0,05,  whFch  yFeFds  a sFgnaF-to-iioise  ratFo  of  20  Foo.  ■ 26  dB 

"lO  iO!) 

The  resuFt  of  the  Feast  squares  Fnverse  fFlter  Fs  shown  Fn  FFg.  h.  The 
spectraF  densFty  of  the  sFgnal  S^(u)  was  assumed  to  l>e  a raFsed  cosFne. 

The  fFFter  Fs  not  abFe  to  resoFve  the  orFgFnal  two  pea  Fes. 

The  resuFt  of  the  project F on  algorFthm,  when  no  restrFctFon  Fs  pFaced 
on  the  sFgnaF,  Fs  gFven  Fn  FFg.  5a  and  5b  for  FO  and  '-O  FteratFons,  respec- 
tFvely.  The  pealcs  are  resoFved;  however,  there  Fs  consFderable  error  outsFde 
the  FntervaF  cor.taFnFng  the  peaies.  When  the  sFgnal  after  each  FteratFon  Fs 
restrFcted  to  onFy  posFtFve  vaFues,  f * (x)  >0,  the  resolutFon  of  the  two 
pealcs  Fs  greatly  Improved.  And  the  overall  error  Fs  reduced  considerably. 

The  results  after  FO  and  kO  Iterations  Fs  given  Fn  FFg.  6a, b. 

The  projection  algorFthm  Fs  versatile  In  that  Ft  Fs  capable  of  utilizing 
a priori  Information  about  the  signal.  If  the  values  outsFde  the  region  of 
the  two  rectangular  pulses  are  known  to  be  zero,  the  algorFthm  can  use  this 
Information.  The  result  of  setting  f(x)  - 0 after  each  subFteratlon  for 
points  outsFde  the  FntervaF  Fs  shown  Fn  FFg.  7a,  b for  1 and  kO  FteratFons. 

The  first  FteratFon  Fs  shown  because  the  algorithm  was  able  to  produce  two 
peaks,  and  this  was  not  the  case  before.  The  fourtleth  FteratFon  Fs  actually 
not  significantly  different  than  the  tenth  Item. on,  which  Fs  not  shown.  A 
more  proper  way  to  utilize  a priori  Information  would  be  to  Incorporate  Ft 

Into  the  structure  of  the  algorFthm.  This  approach  will  be  considered  Fn  the 
near  future. 

The  computation  time  required  for  the  Fnverse  filter  Fs  about  0.8  sec- 
onds. One  FteratFon  of  the  projection  algorFthm  takes  approximately  0.4  seconds. 


The  results  for  the  two-dimensional  Images  indicate  that  the  projection 
method  can  improve  on  the  degraded  image,  but  the  effects  of  noise  dominate 
the  image  after  a certain  number  of  Iterations,  The  optimum  number  of 
iterations  Is  a subjective  quantity.  The  computation  time  required  for  one 
Iteration  on  a 128x128  image  array  is  about  12  seconds. 

iV.  CONCLUSIONS 

The  results  show  that  the  projection  algorithm  can  be  a useful  Image 
restoration  technique.  The  method  is  capable  of  Incorporating  restrictions 
on  the  iignal  and  a priori  knowledge.  The  method  has  been  applied  success- 
fully to  one  and  two-dimensional  signals. 

Future  work  in  this  area  will  involve  the  application  of  the  algorithm 
to  Images  with  more  Involved  types  of  degradation.  The  effect  of  error  In 
the  coefficients  of  the  degrading  transformation  will  be  Investigated. 
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Fig.  ^ Restoration  by  loast'square  Inverse  filter 
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Fit.  7a  Iterative  restoration,  shape  constraint 
(1  Iteration) 
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Fig.  7b  Iterative  restoration,  shape  constraint 
(1*0  Iterations) 
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( SATELLITE  IMAGE  NOISE  REMOVAL 

I 

! 0.  't.  Mitchell  and  P.  L.  Chen 

1 i 

j INTRODUCTION 

I We  are  using  homomorphic  filtering  techniques  to  remove  multiplicative 

noise  effects  such  as  cloud  and  atmospheric  disturbances  In  ERTS  Imagery^  A 
simulated  Implementation  using  simple  computer  generated  pictures  has  been 
1 performed.  The  picture  shown  In  Figure  1 was  multiplied  by  various  types  of 

\ noise  and  then  filtered  to  recover  the  original  picture. 

RANDOM  AND  LOW-FREQUENCY  NOISE 

Low  frequency  noise  was  generated  by  passing  computer  generated  random 
noise  through  a two  dimensional  low-pass  filter.  The  original  noise  was 
uncorrelated  (white)  and  uniformly  distributed  between  0 and  I.  The  picture 
size  Is  6kx6k  pixels.  An  Ideal  rectangular  low-pass  filter  was  used  which 

I 

passed  the  first  8 frequency  components  In  both  the  horizontal  and  vertical 
directions  and  rejected  the  upper  24  frequency  conponents  In  each  direction. 

Figures  2(a)  and  2(b)  show  the  histograms  of  the  random  noise  and  the 
low-pass  noise,  respectively.  The  Intensity  picture  and  the  autocorrelation 
of  the  low-pass  noise  are  shown  In  Figures  2(c)  and  2(d),  respectively. 

HOMOMORPHIC  FILTERING  PROCESS 

Tne  two  dimensional  filtering  to  reduce  the  noise  effect  Is  accomplished 

by  weighting  each  point  In  the  Discrete  Fourier  Transform  of  the  logarithm 

of  the  picture  Intensities  by  a noncausal  filter  function  of  the  form 
Si 

where  Sj  and  Nj  are  the  values  of  the  signal  and  noise  power  spectrums 
at  the  particular  frequency  point  being  considered.  Following  the  filtering 
process  we  take  the  Inverse  DFT  and  then  exponentiate  to  obtain  the  filtered 
picture. 
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Figure  3(a)  shows  the  original  picture  with  mul tipl  Ici^tlve  low-pass 

I 

noise  and  Figure  3(b)  shows  the  resi ’ting  filtered  picture.  Figure  3(c)  shows 
I the  original  picture  with  multiplicative  white  noise  and  Figure  3(d)  Is  the 

resulting  filtered  picture. 


( 


Fig.  k Block  diagram  of  homomorphic  image  filtering 
We  are  now  extending  this  model  to  ERTS  Imagery.  The  power  spectrum  of 
the  signal  and  noise  must  be  approximated  since  they  are  not  directly  avail- 
able. The  signal  spectrum  will  be  approximated  by  measuring  clear  ERTS  frames; 
the  noise  spectrum  wl 1 1 be  obtained  using  multlspectral  channel  classification 
^ on  the  noisy  picture  to  produce  a binary  picture  of  the  noise. 
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INTRODUCTION 


A Kalman  filter  for  restoration  of  Images  degraded  by  any  LSI  (linear 
shift  Invariant)  point  spread  functions  Is  developed  here.  Also  mentioned 
here  are  certain  flaws  In  some  earlier  work  done  by  other  authors  [1,2]  In 
developing  Kalman  filtering  techniques  for  Image  enhancement. 


II.  IMAGE  HODEL 


The  Image  is  assumed  to  be  a two-dimensional  wide-sense  Markov  field 
[3]  with  exponential  autocorrelation  function,  and  defined  on  a discrete 
domain  as 


s(i,j)  - p s(I-l,j)  + p.s(l-l,j-l)  + p.s(l,j-l)  + oj(l,j) 


0) 


whe  re , 


s(l,j)  “ the  undlstorted  Image  at  the  point  (l,j) 


Py  t ■ vertical  and  horizontal  correlation  coefficient 


"d  ■ ■'’yPh 


* A sequence  of  zero  mean  uncorrelated  random  variables. 


The  degraded  observation  y^(I,j)  Is  given  by 

I 


:i 


kf 

^ 9'(K'A')  + n'(l,J)  (2) 

" •-“o  * 


where 


1 


I 


1 


q'{y'  ,i-')  ■ the  point  spread  function,  assumed  zero  outside  of 
-kQ^  K'<k'  and  "tj  1 < Af 

n'*(l,J)  - A sequence  of  zero  mean  uncorrelated  noise  with  variance 

o2 

n 


We  shall  now  define 


y(i,j)  - y'(i-kJ.J-Ao) 


P - k"  + kj 


q - r + i; 


k ■ + k^ 

K T Nq 

I m I'  + i.' 

n(l,J)  - n"(l-kj  , J-ftJ) 


g(p,q)  • g‘'(k'',J.‘‘)  - g‘*(p"kj  , q-^p 


then  we  get  a shifted  version  of  the  observation 

k A • 

y(l,j)  - E I s(l-p,j-q)  g(p,q)  +n(l,j) 

p-0  q«0 


(3) 


1 


and. 


E[n(l  ,j)n(k,^')]-  A(l-k)A(j-i-) 


where,  A ■ Kronecker  delta. 

If  the  size  of  the  observed  Image  y(l,j!  Is  mxn  then  we  let 

s(l)  - [s(l,D  s(l,?.)  ....  s(l,n)]'’’ 
y(i)  - [y(i,D  yd. 2)  — y(i,n)]^ 
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w(I)  - [(ii(I+l,l)  wd+1.2)  ....  wO+l,n)]T 


The  n we  have 


X(l+1)  - A X (i)  + B 0)(I) 

y(i)  - C X (i)  + n(i)  W 

C ■ 1 (k-0  arid  2.-0)  give*,  the  observation  corrupted  by  additive  noise  only. 
III.  MOTIVATION  FOR  THE  MODEL 

As  will  be  seen  In  later  sections  the  model  described  In  eqn.  {k)  will 
lead  to  a recursive  restoration  technique  In  which  the  estimate  of  the  Image 
s(I,j)  Is  expressed  I.t  terms  of  past  estimates  and  the  obse'"vatlons  along  a 
complete  row  of  the  picture.  Hablbl  and  Nahl  [1,2]  have  ured  the  model 
described  by  eqn.  (1)  In  enhancement  of  Images  corrupted  by  additive  noise 
only,  and  have  developed  an  estimate  s(l,j)  which  is  expressed  In  terms  of 
previous  estimates  and  the  observation  y(l,j)  only  at  the  point  (l,j).  But, 

It  can  he  shown  that.  In  such  a formulation  the  estimation  error  at  the  point 
e(i*j),-  Is  not  necessarily,  orthogonal  to  the  observation  y(p,q)  for 
l£  p<  1-1,  q-j  and  for  p-1,  1 q ^ j"<.  In  this  sense  the  estimate  Is  sub- 
optimal  though  computationally  fast  since  it  uses  the  observation  only  at  one 
point,  rather  than  on  a complete  row.  It  can  also  be  shown  that  when  this  is 
extended  to  estimation  of  images  degraded  by  LSI  point  spread  function  the 
suboptimality  worsens.  The  recursive  estimates  are  functions  of  estimates  In 
the  "past"  and  the  observations  at  the  "present".  It  was  suspected  that  the 
above  mentioned  problem  may  be  due  to  the  fact  that  In  a two-dimensional 
Markov  field  the  "present"  (as  oppos'id  to  "past"  and  "future")  Is  a curve 
[see  e.g.  k,S]  rather  than  just  one  point  as  In  one-dImens lonal  Markov 
sequences.  The  model  described  by  eqn.  {k)  Implicitly  defines  the  "present" 
as  one  cfxnplete  row  rather  than  just  one  point  and  indeed  this  eliminates  the 
suboptimality  mentioned  above. 
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iv;  THE  FILTER/PREDICTOR 


We  now  can  derive  the  recursive  MMSE  estimate  of  X(I)  In  terms  of  esti" 
mates  X(P),  I and  y(I)  by  following  siml  lar  d«rTvatlon  for  the 

estimate  of  a continuous  time  domain  process  as  In  [6].  If  Instead  wewanted 
to  estimate  X(I)  In  terms  of  X(P),  I <P  and  y(l-l)  we  could  use  the 

resuK  directly  from  [?].  Both  the  results  are  given  below. 

X{U\)  - A x"(I)  + H(I+l)[y(I+I)  - CAX(l)] 

H(i+i)  - Q(i)  c^[cQ(i)c^  + 

Q(i)  - AP(i)A^  + BZjo 

. Q(I+I)  - A[l  - H(I  + I)C]  Q(1)A^  + BEiob''^  (5) 

And 

X(I+I)  - AX(i)  + H(I)[y(l)  - CX(1)] 

H(I)  - AP(I)c''^[CP(I)c''^  + 

P(I+1)  - [A  - H(I)C]  P(i)[A  - H(I)C]'^  + BZo)  b"^  + H(I)zy(l)  (6) 

where,  - 

L -E[u)(I)J(i)l  . - Eln(I)n^(I)]  , P(I)  - E[e(I)e' (l)] 

n n 

Compare  the  result  given  In  eqn.  (5)  with  the  corresponding  result  reported 
In  [8]. 

A 

The  estimate  of  5(i,J)  Is  obtained  by  taking  the  jth  row  of  X(I).  It  Is 
necessary  to  have  the  Initial  condition  P.{I)  to  determine  H(I),  which  Is 
analogous  to  the  coefficient  n(l,j)  In  the  earlier  work  [9l  reported  by  these 
authors.  It  Is  expected  that  just  as  with  ri(I,j)  In  [9l  the  effect  of  Initial 
condition  on  H(i)  will  diminish  as  I Increases.  Also,  because  of  modifica- 
tion In  the  definition  of  "present"  the  smearing  of  edges  [see  9l  Is  expected 
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to  be  less  severe.  In  the  presence  of  only  LSI  degradation  and  no  additive 
noise  H(i)  anl  X(i)  In  eqn.  (5)  reduce  to 


H(I) 


X(I)  - H(i)y(I)  (8) 

This  gives  the  inverse  filtering  technique  without  going  to  the  frequency 
domain  and  having  to  deal  with  zeroes  In  the  transform  of  the  point  spread 
function.  The  sensitivity  of  the  inverse  filtering  to  any  holse,  however, 

Is  expected  to  be  high. 

It  must  be  pointeo  jut  that  the  matrices  C,  A,  (i,  etc.  can  be  very  large 
even  for  a moderate  size  picture  (e.g.  128x128)  and  point  spread  function 
(e.g.  10x10).  In  such  a case  the  picture  may  be  brcken  to  several  smaller 
size  pictures  concatenated  and  process  each  of  the  smaller  pictures  separately 
and  concatenate  the  pieces  back  together.  Similar  work  Is  being  done  for 

linear  shift  variant  degradation  of  Images. 

V.  THE  PROOF  OF  EQN.  (5) 

Based  on  the  model  described  by  eqn.  (4)  we  want  the  linear  HHSE  estimate 

A 

X(I)  of  X(i)  in  terms  of  y(l) y(i),  i.e.. 


X(i)  - Z H(i ,p)y(p) 

p-1 


wl  th 


H(i ,p)  such  that 


E[{X(i)  - X(i)}y'’’(q)]  - 0 q < I 

T ^ T 

E[X(l)y  (q)]  - Z H(i ,p)E[y(p)y  (q)] 

p-i 

E[X(I+l)y^(q)  - E[(AX(I)  + Bw(l)}y^(q)] 
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- AE[X(l)y^(q)  + B E[w(l)y^(q) 


Since,  y(q)  ■ CX(q)  + n(q) 

- CAX(q-l)  + CB(jj(q-l)  + n(q) 

y(q).  aJ*  q ^ contain  oj(l)  and  hence  from  eqn.  (1) 

E[a)(l)y^(q)l  - 0 

E[X(i+l)y'^(q)]  - S AH(l,p)  E[y(p)y'^(q)l 
P-1 

But  from  eqn.  (10) , 

_ 1 ♦•1  T 

E(X(l+l)y’(q)  - I H(1+1,P)  E(y(p)y  (q)l 
P“1 

- Z H(l+l,p)E[y(p)y^(q)j+  H(l+l,l+l)E[y(l-l)y^(q)  ] 

P-1 

E(y(l+1)y^(q)l  - E[{CX(1+1)  + n(l+1))y^(q)] 

- CE[X(i+l)y^(q)  + E[n(l+l)y^(q)l 

Since,  for  all  q £ 1»  y(q)  does  not  contain  n(l+l)  E[n(l+1)y  (q)l  - 0. 

The re fore, 

T ^ T 

E[y(l+l)y  (q)l  - ^ CAH(l,p)  E[y(p)y  (q)] 

P“1 

E(x(l+l)y^(q)]  - ^ [H(l+l,p)  + H(l-H,l+l)CAH(l,p)E[y(p)y^(q)]  (12) 
P-1 

Combining  eqn.  (11)  and  (12) 

I [(A  - H(l+l,l+l)CA)H(l,p)  - H(l+l,p)]E[y(p)y^(q)l  - 0 
P-1 


or,  ZH^(I,p)  Ety(p)y^(q)l  - 0 

P-1 

where  H^(I,p)  - (A  - H(l+1 ,1+l)CA)H(l,p)  - H(I+l,p) 


Let  X"(I)  - E {H(l,p)  +H"(l,p>}y(p) 

p-1 


E[{X(i)  - r(I))y^(q)]  - E[X(l)y^(q)]  - Z H(l  ,p)E[y(p)y”^(q)] 


SH'(I,p)E[y  (p)y  (q)  ] 

p-1 


- 0 from  eqn.  QO)  and  (15) 


. . Both  X(l)  and  X‘'(l)  are  linear  HHSE  estimates  of  X(l).  Hence,  from 
uniqueness  of  the  linear  HHSE  estimate 


A A 


X(l)  - X"(l) 


A A 


or  X(l)  - X"(l)  - Z H"(l,p)y(p)  - 0 

P-1 

Since  eqn.  (14)  Is  true  for  all  values  of  Y(P)  It  roust  be  true  that 
H"(l,p)  - (A  - H(l+l,l+l)CA)H(l,p)  - H(l+l,p)  - 0 


X(l  + 1)  - Z H(l  + l,p)y(p) 

p-1 

i 

- Z (A  - H(l+l,l+l)CA)H(l,p)y(p)  + H(l+1 , 1+1 )y (1+1 ) 

p-1 


- (I  - H(I+1,I+1)C)AX(1)  + H(l+l,l+l)y(l+l) 


t 


I"  ~ 


fi  ' 

- AX(i)  + H(i+l)[y(l  + l)  - CAXO)]  (l6) 

by  writing  H(i  + l)  for  H(I+1,I-H). 

Hence,  the  error  e(i+l)  - X(I+1)  - X(l+1) 

- AX(i)  + Bo)(I)  - [AX(l)H(l+l){CAX(i)  + CBo)(i)  + n(I+l)  - CAX(i)}] 

- (I  - H(l+l)CAe(I)  + (I  - H(l+l)C)Bo)(i)  H(i+l)n(l+l)  (17) 

Since,  E[e(l+1  )y^(q)  ] 0 for  all  q 1+1 

E[e(I+l)  r (I)]  - 0 
.*.  E[e(I+l{y(l+l)  - CAX(I)}^]  - 0 

or  E[e(I+l){CAe(l)  + CBu>(i)  + n(l+l)}"*^]  ■ 0 (18) 

Since  E[w(I)y^(q)]  ■ 0 q £ • 

E[u(I)e^(I)]  - 0 

Substituting  (17)  In  (I8)  and  writing  E[e(l)e^(l)]  - P(l)  we  get, 

(I  - H(i+l)C)AP(i)A^C^  + (I  - H(I+1)C)BE  bV  + H(l+l)Z  - 0 

a>  n 

or  H(I+1)  - (AP(I)A^  + BE  B^)C^[C(A  P (I )A^  + BE  B^)C^  + E ]’^ 

Oi  0)  n 

- QC''’[CQC''’  + E^]’‘  (19) 

P(I  + 1)  - E[e(I  + I)e''’(I+l)] 

- E[e(l+l)x''’(l+l)] 

- E[e(l+1){AX(I)  + Bo)(l)}'*’l 
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Substituting  eqn.  (17)  we  have 


P(I+1)  - ( 1 - H(I+1)C}  APCDa'*^  + { I - H(I+1)C}  BEy 
- [I  - H(l+l)Cl  d(l) 
or  Q(I+1)  - APd+OA*^  > BL^  b”^ 

- All  - H(1+1)C](1(1)A^  + BE^  b"^ 

Putting  together  eqn.  (16)»  (I9)i  and  (20)  we  get  eqn.  (5)* 
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CHARACTERIZING  THE  INFORMATION  CONTENT  OF  MULTISPECTRAL  IMAGES 

p.  A.  WIntz 

I.  INTRODUCTION 

The  information  contained  in  mul t i spectral  imagery  of  the  earth's  sur- 
face is  spread  through  four  dimensions  - 2 spatial,  1 spictral,  and  1 tem- 
poral. Whether  the  information  is  extracted  by  a person  or  a machine  the 
magnitude  of  the  information  extraction  task  could  be  lessened  if  the  data 
bulk  could  first  be  reduced.  In  the  case  of  machine  processing  there  are 
three  sources  of  savings  - the  time  required  to  Input  the  data  into  the 
machine,  the  amount  of  memory  required  to  store  the  data,  and  the  amount  of 

machine  time  required  to  extract  the  information. 

The  problem  is  to  minimize  the  cost  of  extracting  Information  from  the 
data  under  the  constraint  that  the  degradation  in  the  performance  of  the 
information  extractor,  e.g.,  classification  accuracy,  not  exceed  a prescribed 

amount. 

Some  source  coding  algorithms  may  be  reasonably  efficient  for  our  pur- 
pose. An  algorithm  that  appears  to  serve  both  purposes  in  some  applications 
is  the  principal  components  (eigenvector)  transformation.  Fig.  1 shows  the 
results  on  an  experiment  performed  on  some  12  channel  aircraft  data.  The 
right  hand  scale  gives  the  mean  square  error  between  the  original  and  re- 
constructed images  and  the  left  hand  scale  lists  the  percent  classification 
accuracy  obtained  by  inputing  the  transformation  coefficients  directly  into 
the  pattern  recognition  machine.  Fig.  la  resulted  from  using  the  transform- 
ing 10x10  sub  arrays  in  the  spatial  domain.  Compression  ratios  of  4 in  the 
spectral  domain  and  10  in  the  spatial  domain  were  achieved  before  classifi- 
cation accuracy  degraded.  Since  the  number  of  computations  required  to 
implement  the  maximum  likelihood  algorithm  for  Gaussian  data  Is  proportional 
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to  the  square  of  the  dimensionality  of  the  data,  a significant  reduction  in 
the  computational  load  of  the  pattern  recognition  machine  was  achieved. 

Any  solution  to  this  problem  must  be  both  data  and  data  user  dependent. 

The  data  dependence  Is  similar  to  the  data  dependence  of  source  codes:  they 

must  be  matched  to  the  characteristics,  structure,  statistics,  etc.  of  the 
data.  The  dependence  on  the  data  user  and  his  data  requlremeflts  can  be 
illustrated  by  looking  at  two  extreme  cases:  If  the  user  data  requirements 

are  not  well  defined  (we  do  not  know  a priori  who  the  data  users  will  be  or 
what  they  will  require  of  the  data)  then  the  only  satisfactory  solution  may 
be  to  ailow  no  distortion  of  the  data,  i.e.,  use  an  information  preserving 
code.  At  the  other  end  of  the  scale  suppose  that  we  know  that  all  potential 
users  will  do  spectral  signature  analysis  and  that  we  also  know  the  classifi- 
cations. Then  we  could  classify  each  resolution  element  and  code  the 
classification  map.  If  K classes  are  required  then  log^K  bits  per  resolution 
element  are  needed  to  code  the  classification  map.  If  the  original  data 
consisted  of  N multispectral  channel  and  M bits  per  resolution  element  then 
the  original  data  set  contained  NM  bits  per  resolution  element.  For  example 
EkTS  data  has  NM=28  whereas  a classification  map  with  32  classes  would 
require  log^  32=5  bits  per  resolution  element.  The  classification  map  could 
itself  be  coded  by  contour  coding,  the  Duan-WIntz  [A]  folded  DPCM  method, 
etc.,  to  achieve  fewer  bits  but  the  cost  of  decoding  would  probably  be  pro- 

hi bi tive. 

Consider  an  intermediate  case.  Suppose  that  the  set  of  users  are  con- 
strained to  spectral  signature  analysis,  but  the  classification  categories  are 
unknown.  Can  we  provide  the  data  user(s)  with  a significantly  reduced  set  of 
data  in  a format  that  minimizes  the  amount  of  computation  he  must  do  while 
allowing  arbitrary  classification  categories  with  negligible  degradation  in 
classification  accuracy? 


17<4 


F , 

E 

% 

r 

1 

-•  i. 

1 

J 

i-  ■ 

S.  1 

f.  F 

I 

r 

i-  i 

? 

r 

} 

r 

f j 

' 

(c)  Channel  3 


(d)  Channel  ^ 


Fig.  2 Sections  of  a frame  of  ERTS-I  MSS  Data 
(LARS  Run  No.  72032806) 
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II.  ADAPTIVE  FEATURE  SPACE  CODING 


Our  approach  to  this  problem  Is  to  determine  all  the  distinguishable 
clusters  In  the  spectral  dimension  of  the  data  and  to  code  the  clusters.  That 
Is,  we  first  use  a clustering  algorithm  to  determine  the  clusters;  then  we 
determine  the  mode  center  of  each  cluster;  then  we  label  each  mode  center; 
finally,  we  generate  a cluster  map,  l.e.,  a map  with  each  resolution  element 
labelled  with  Its  cluster  number.  Since  every  conceivable  category  cor- 
responds to  one  of  the  clusters,  any  data  user  can  classify  the  data  Into  any 
Set  of  categories  using  a table  look  up  procedure.  That  is,  he  first  deter- 
mines which  set  of  clusters  correspond  to  which  category  (training)  to  generate 
the  table,  and  then  he  classifies  each  resolution  element  by  using  Its  cluster 
label  to  address  the  table.  For  the  ERTS  frame  presented  In  Fig.  2 we  found 
96  distinguishable  spectral  domain  clusters  so  that  the  cluster  map  required 
7 bits  per  resolution  element  whereas  the  original  ERTS  data  required  28  bits 
per  resolution  element.  Every  conceivable  classification  category  corresponds 
to  one  or  more  of  the  clusters  so  that  any  data  user  can  classify  any  resolu- 
tion element  by  mapping  its  cluster  label  into  the  appropriate  category  label. 

As  an  example,  we  used  this  method  on  data  from  channels  2 and  3 of 
the  ERTS  frame  presented  in  Fig.  2.  With  a window  of  6x6  resolution  elements 
around  each  mode  center  In  the  2-dlmensional  feature  space  our  algorithm  found 
96  distinguishable  clusters  so  that  the  cluster  map  required  7 bits  per 
resolution  element.  Applying  the  Duan-Wintz  folded  DPCM  algorithm  to  the 
cluster  map  reduces  this  to  ^ bits  per  resolution  element.  With  this  cluster 
map  as  the  Input  data  to  the  classification  algorithm  with  the  three  classes 
corn,  soybeans,  and  others  we  achieved  a classification  accuracy  of  11%.  The 
standard  LARSYS  maximum  likelihood  algorithm  operating  on  the  original  data 
achieved  a classification  accuracy  of  78^  for  the  same  classes. 
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Let  us  now  look  at  this  method  and,  in  particular,  the  example  of  the 
preceding  paragraph  ^rom  a different  point  of  view.  The  2-dlmens lonal 
feature  space  obtained  by  plotting  the  channel  2 Intensity  against  the 
channel  3 Intensity  consists  of  127x127  points  and  each  ERTS  resolution  ele- 
ment falls  into  one  of  these  points.  Our  method  achieves  a data  bulk  re- 
duction by  requantizing  this  space  more  coarsely.  This  is  accomplished  by 
first  making  a scatter  plot  of  the  data,  then  moving  a 6x6  window  around  to 
Mnd  the  location  having  the  largest  number  of  data  points  within  the  window. 

This  6x6  subspace  of  the  feature  space  is  then  labelled  ff\  and  all  data  points 

within  the  window  are  given  this  label.  This  corresponds  to  constructing  a 
2-dimensionai  6x6  quantization  bit  at  this  location  and  rounding  off  all  data 
within  the  bin  to  a common  quantization  level.  Next,  the  6x6  window  again 
moved  to  again  find  the  location  having  the  largest  number  of  data  points 
within  the  window  only  this  time  with  the  constraint  thac  It  cannot  overlap 
quantization  bln  §2.  This  seco: d 6x6  array  is  labelled  #2.  This  process  Is 

repeated  until  we  run  out  of  data  points.  For  the  ERTS  frame  we  ended  up  with 

96  bins  ordered  according  to  the  number  of  data  points  within  each  bln.  Some 
regions  of  the  feature  space  contained  no  data  points  and  so  no  bins  were 
constructed  in  these  regions.  In  general,  there  are  some  gaps  between  bins 
that  contain  data  points.  Refer  to  Fig.  3 and  suppose  that  we  first  found  bin 
#1,  then  bin  /'2,  then  bin  §2 . Due  to  the  constraint  that  bins  can  not  overlap 
gaps  can  occur  between  the  bins.  Data  points  within  the  gaps  are  assigned  the 
label  of  the  nearest  bin,  i.e.,  rounded  off  to  the  nearest  bin.  This  Is 
equivalent  to  letting  the  bins  grow  to  take  in  these  neighboring  points  so  that, 
finally,  some  of  the  quantization  bins  are  irregular  In  shape  (not  6x6).  For 
the  data  of  Fig.  2 this  2-dimensional  coarse  requantization  of  the  data  led  to 
a 1%  degradation  in  the  classification  accuracy. 


178 


III.  NONADAPTIVE  FEATURE  SPACE  CODING 


Another  way  to  more  coarsely  quantize  the  2-d Imens ional  feature  space  Is 
to  fill  the  feature  space  with  a regular  array  of  6x6  quantization  bins  as 
illustrated  in  Fig,  k.  This  corresponds  to  more  coars  *y  quantizing  both  of 
the  channel  intensities  and  would  correspond  to  returning  few^r  significant 
bits  for  each  for  the  natural  code  If  the  array  size  were  a power  of  two. 
Whereas,  the  method  of  Fig.  3 is  data  dependent  (adaptive)  the  method  of 
Fig.  k is  not.  Obviously,  the  method  of  Fig.  ^ is  significantly  easier  to 
implement.  When  this  method  was  used  on  the  ERTS  frame  of  Fig.  2 a classifi- 
cation accuracy  of  73?'  was  achieved. 

IV.  STAGGERED  DPCM 

One  method  for  reducing  the  number  of  bits  required  to  characterize  a 
mul tispectral  image  is  to  take  advantage  of  the  spatial  and  spectral  redun- 
dancy by  simply  neglecting  some  of  the  pels.  In  the  reconstruction  process 
the  missing  picture  elements  are  estimated  from  the  retained  picture  elements 
and  inserted  into  the  reconstructed  image..  The  retained  picture  elements  are 
coded  by  DPCM.  To  minimize  the  loss  in  spatial  and  spectral  resolution  the 
retained  picture  elements  from  successive  lines  are  staggered,  e.g.,  if  the 
even  numbered  pels  are  deleted  r *om  the  first  row,  then  the  odd  numbered  pels 
are  deleted  from  the  second  row,  etc.  so  that  each  deleted  picture  element  Is 
surrounded  on  ail  four  sides  by  retained  picture  elements  and,  in  the  re- 
construction process,  can  be  estimated  from  its  four  nearest  neighbors.  The 
deleted  pels  can  also  be  staggered  in  the  spectral  dimensions  as  illustrated 
in  Fig.  5 so  that  half  the  data  points  in  each  spectral  vector  are  retained.' 

A range  of  data  compression  ratios  can  be  obtained  by  varying  the  number  of 
bits,  i.e.,  the  number  of  quantization  levels,  used  to  code  the  pel  differ- 
ences. 


ISO 


Fig.  5 


This  technique  was  tested  on  two  ERTS  frames,  LARS  run  no.  730^1801  pre" 
sented  In  Fig.  2 and  LARS  run  no.  72032806  presented  In  Fig.  6.  The  recon- 
structed data  was  Inputed  to  the  standard  LARSYS  maximum  likelihood  classifi- 
cation algorithm  and  the  resulting  classification  accuracies  are  presented  In 
Fig.  7 for  a range  of  compression  ratios.  The  degradation  In  classification 
accuracy  relative  to  the  accuracy  obtained  from  the  original  data  (which 
corresponds  to  a compression  ratio  of  I In  Fig.  7)  is  surprisingly  small  out 
to  compression  ratios  of  7. 

We  also  generalized  the  same  basic  strategy  to  retain  every  third  pel  and 
the  performance  of  this  strategy  Is  also  plotted  In  Fig.  7 at  a compression 
ratio  of  9.  It  appears  that  this  does  not  work  as  well  as  retaining  every 
second  pel  and  using  a coarser  quantizer  to  achieve  the  same  compress  Ion  rat  lo 
Also  plotted  on  Fig.  7 are  the  performances  of  the  two  ciuster  coding 
techniques  discussed  In  the  preceding  sections. 
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Sections  of  a frame  of  ERTS  - I MSS  Data 
tcARS  Hun  No,  730^^801) 


J 
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We  fInH  It  surprising  tf *t  the  simple  data  deletion  technique  performs 
so  well.  It  Is  not  likely  that  a slmplor  method  could  be  found  that  performs 
this  well  and  so  this  appears  to  be  an  optimum  solution  to  the  problem  of 
Image  coding  relative  to  u criterion  based  on  the  classification  accuracy 
achieved  with  the  coded  data.  On  the  other  hand.  It  has  no  utility  to  thj*. 
problem  we  posed  at  the  beginning  since  the  compressed  data  must  be  decoded 
before  It  can  be  classified. 
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Present  local  facilities  available  to  us  include: 

I  POP  11/45  computer;  32K  core;  disk  operating  system  software 

1 20-mi  1 1 Ion  word  disk 

2 9"track  tape  drives;  800  and  1600  bpi  densities 
I /"track  tape  drive;  556  and  800  bpI  densities 

3 CRT  terminals 
I card  reader 

I paper  tape  reader 

I line  printer;  600  lines  per  minute.  Impact  type 
I electrostatic  graphics  hard  copy  system,  200  points  per  Inch 
I video  display  with  disk  storage;  3 color  refreshed  display,  512x512 
points,  8 bl ts/color/point;  graphics  overlay  and  Interactive  capabili- 
ties 

I flying  spot  scanner 

The  hardware  interface  between  our  computer  and  the  IMP  has  been  completed. 
We  are  currently  developing  an  ELF  software  system  for  use  o«  our  system.  We 
are  also  investigating  the  use  of  UNIX  as  a replacement  for  both  DOS  and  ELF, 
Basic  network  data  transfer  capabilities  will  be  available  in  June  1975 
and  complete  hardware  and  software  capabilities  will  be  available  In  late  I975* 
In  addition  to  our  own  facilities  we  have  the  use  of  the  Purdue  University 
Computing  Center's  CDC  6500  and  the  IBM  360/6/  at  the  Laboratory  for  Application 
of  Remote  Sensing.  Batch  picture  processing  programs  with  graphic  output  have 
been  written  for  the  CDC  6500,  and  the  LARSYS  operating  system  at  LARS  allows 
us  to  do  interactive  processing  of  multlspectral  data.  Local  intera'C6*lKMr‘ 
processing  Is  also  available  on  our  POP  11/45  computer.  A software  package 
(AMIDS)  developed  at  RADC  for  pattern  analysis  is  being  transferred  to  our  PDP. 
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MISSION 

'Ranie  Air  Devetopn'isnt  Cetiter 


RADC  iMtIm  principal  KFSC  oxgmniz»tion  chMtg»d  with 
punning  und  ^x^coting  the  USAF  exploretorg  end  edv^d 
development  progreme  for  informetion  acienoee,  intelli 
g^e^Lmnd,  control  end  coemnications 
product*  and  service*  oriented  to  the  need*  of  USAF. 
Primary  MX  mission  areas  are  communications , elect^ 
magnetic  guidance  and  control,  surveillance  of  gnuM 
aiS\eros^ce  objects.  Intelligence  data 
handling,  information  system 
reliability,  maintainability  and 

ha*  mission  responsibility  as  assicried  by  for  de 

monatration  and  acquisition  of  selected  subsystems  and 
systess!  in  the  intelligence,  mapping,  charting,  caemmna, 
control  and  conwurlcwtiotts  area*. 


